{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Discussion of the GP marginal likelihood upper bound\n",
    "--\n",
    "\n",
    "*Mark van der Wilk, August 2017*  \n",
    "*Small edits by James Hensman 2017*\n",
    "\n",
    "See [gp_upper](https://github.com/markvdw/gp_upper) for code to tighten the upper bound through optimisation, and a more comprehensive discussion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T12:31:18.299482Z",
     "start_time": "2018-06-19T12:31:17.188860Z"
    },
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "%matplotlib inline\n",
    "matplotlib.rcParams['figure.figsize'] = (12, 6)\n",
    "plt = matplotlib.pyplot\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import gpflow\n",
    "from gpflow.test_util import notebook_niter\n",
    "\n",
    "import logging\n",
    "logging.disable(logging.WARNING)\n",
    "\n",
    "from FITCvsVFE import getTrainingTestData"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T12:31:18.303116Z",
     "start_time": "2018-06-19T12:31:18.300562Z"
    },
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "X, Y, Xt, Yt = getTrainingTestData()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T12:31:18.316919Z",
     "start_time": "2018-06-19T12:31:18.304403Z"
    },
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def plot_model(m, name=\"\"):\n",
    "    pX = np.linspace(-3, 9, 100)[:, None]\n",
    "    pY, pYv = m.predict_y(pX)\n",
    "    plt.plot(X, Y, 'x')\n",
    "    plt.plot(pX, pY)\n",
    "    try:\n",
    "        plt.plot(m.feature.Z.value, m.feature.Z.value * 0, 'o')\n",
    "    except AttributeError:\n",
    "        pass\n",
    "    two_sigma = (2.0 * pYv ** 0.5)[:, 0]\n",
    "    plt.fill_between(pX[:, 0], pY[:, 0] - two_sigma, pY[:, 0] + two_sigma, alpha=0.15)\n",
    "    lml = m.compute_log_likelihood()\n",
    "    plt.title(\"%s (lml = %f)\" % (name, lml))\n",
    "    return lml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Full model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T12:31:18.930900Z",
     "start_time": "2018-06-19T12:31:18.318308Z"
    },
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEICAYAAABLdt/UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJztnXl83HWd/5+fuWdyt0mbkqMXpbSUlkAtlauCoChQQESxCgLranVZ111dFfnpanc9Vndd2fVgXXdxV6gIAnLIVRBUjraUhpbe0CtNeiVt7smc38/vj5lJZzIzyWSOzJH38/HIo+n3/Hwz831939/3530orTWCIAhC6WDK9wAEQRCE7CLCLgiCUGKIsAuCIJQYIuyCIAglhgi7IAhCiSHCLgiCUGKIsAtpo5T6tVLquvDvtyqlXk7zOO9RSrVnd3TCSJRSdqXULqVUXb7HIuQWEXYhLZRSi4ElwGP5HkuqKKX+RSn1tlKqPyxwt0Stq1VKvaKUOqGU6lFKvaaUunCUYzUopR5TSp1USrUrpVaPWG9WSv2TUupw+HytSqnqqPVzlFJPhtd1KaW+P2L/m5RSO5VSg0qpvUqpi8PLZymltFJqIOrn61H7TVFK/SZ8HV1KqfuVUpUAWmsv8D/AVzP9WwqFjQi7kC6fAe7XxZXhNghcA1QBnwTuVkpdEF43ANwO1AE1wD8DTyilLEmOdR+wH5gOXAV8Ryl1adT6bwEXAO8GKoGbAQ+AUsoGrAP+ANQDjeHjEV5/Rfj8twEVwCXAvhHnr9Zal4d//jFq+T+Fxz8bmBse3zej1q8FPqmUsie5LqEEEGEX0uUDwB+TrQxblZ+LspD/USk1Vyn1qlKqTyn1YFjgJgyt9T9orXdprQ2t9Qbgz4SEF621R2u9W2ttAAoIEhLIKQmurRx4D/BtrbVfa70F+C2hBwNKqRrgC8Bfaq0P6hDbtNae8CFuBQ5rrX+otR4Mn3tr1Cm+BazRWq8Pj7VDa92R4mXOBn6nte7TWvcCjwJnRf0N2oFuYHmKxxOKEBF2YdwopcoICcjuMTZ9P3AeIRH5MvBz4BNAE7AI+FiK59sado8k+vlpmtfgBN4FbB95LkKW9ePAL7TWxxPtPuLfyO+Lwr+fDQSADyuljiql9iil/ipq2+XAAaXU02F3yUtKqbPD5zcDS4E6pdQ7YTfPj8PjjeZgeN29SqnaqOU/Aa5WStWEHzA3AE+P2HcnITeaUKKIsAvpEPEV94+x3ffDluN2YBvwnNZ6X9iSfBpoSeVkWuvFWuvqJD+fS/Ma7gG2AM+OPBch18kqIOFksNa6H3gF+LpSyqGUOpeQgLrCmzQScvecQegB+GHgm2EXS2T9TcC/A6cBvwceC7/BTAes4X0uBs4h9Hf6f+F9uwg9kGYSemhWAPdHDW8zYANOhH+CwMiHXz+nPkOhBBFhF9KhJ/xvxRjbHYv6fSjB/8uzOaholFL3RE0ufm3Euh8Qsq4/kmiOIOwa+TXwVaVUMsv244RE+xDwM0I+8khkz1D43zVa66Gwm+UB4INR61/WWj+ttfYB/wJMBRZE7fsfWusjWusu4IeRfbXWA1rrTVrrgNb6GHAH8D6lVOSzeBDYQ+izqQT2EuW/D1PBqc9QKEFE2IVxo7UeJCQYZ0zE+ZRS20dEgUT/3JNkjKujJhe/E3WsbxGaH3if1rpvjFNbgTlJjn9Qa3211rpOa30+UAtsDK+O+MujHxrRv28d8f/o43YTekAk2zdul/C/kXv5HOA/w777AUJvJh8csc8CQm8rQokiwi6ky1PAiok4kdb6rCiRHvmzeuwjhFBK3UnIxXK51vrEiHXLlVIXKaVsSimnUuorhNwiG5Ica4FSqiK8/SeA9xGyrNFa7yU0MXtXOHZ8ASHXy5Ph3e8DliulLg/71L9AyMWyM7z+XuCvlVLTwn7yv43sq5Q6Xyk1XyllUkpNJeTOeSns3gJ4HfhU+BqcwKc59aBBKdVAaEJ4fap/N6H4EGEX0uXnwMeVUmrMLQuH7wDNwDsJ3DR2QhOPJ4AOQlbuVVrrwwBKqY8rpaInWt9PKASxG1gNXKm17oxa/zFCfvAThHzoX9davwCgtd5NaBL5nvD+1wIrw24ZgH8kJNB7CIl9K/Dt8Lo5wDOE/OTbAC+xk9C3A7MIWf0d4e0/GbV+FfC/4Zh2oURRxRWGLBQSSqm1wINa69/leyzC2IRj17cAlySJ9hFKBBF2QRCEEkNcMYIgCCWGCLsgCEKJIcIuCIJQYiQrcJRTamtr9axZs/JxakEQhKLljTfe6NJaj1l2OS/CPmvWLDZt2pSPUwuCIBQtSqmDqWwnrhhBEIQSQ4RdEAShxBBhFwRBKDFE2AVBEEoMEXZBEIQSQ4RdEAShxBBhFwRBKDEyFnalVJNS6kWl1I5wQ4S/ycbABEEQhPTIRoJSAPii1npzuD3XG0qpdVrrHVk4tiAIJc6gN8CgL4DWYGiNQmE1K2wWEzaLCbvFnO8hFh0ZC7vW+ghwJPx7v1JqJ9AAiLCngGFoNBApn2wxi3dMmBz0e/z0uP34g8ao25mUwm414bCYsVtDQm82FVN/l4knqyUFlFKzCHVUj2snppT6NKE2XTQ3N2fztAVN0ND4gwZBQxMIavyGQSCoCYT/NUbUwzebQpaKw2KmymnFJF9goQQ5MeCld8if0raG1gz5ggz5gsPLrOaQNW8zm7BaTFhMKvQjhhGQRWFXSpUDDwNfSNQkWGv9c0Lt1Fi6dGna3T36PX6M0R/wWUWP6CMc0eGIla3DyyK/G1oTNDSGAUGtGW8jk6Bx6kvc7wlQU2alwmHNyrUIQiHQO+RPWdST4Q8a+IMGgyOWK6UwK4XJFDKSFAqTCi1XClR4Gwj9Hton8TkU2Teq7FYTDmvuXUtZEXallJWQqN+vtX4kG8dMRiqvbqVCwDDo7Pcy4A1QX+mguNqLCkI8g94AJwZy125Va01AayhQiah22SZE2LMRFaOA/wZ2aq1/mPmQhJEM+YJ05vBmEISJwB80ON4v3+OJIBsOqQuBm4HLlFJvhn8+mIXjClEMeAL0uH1jbygIBUqP2z9u16SQHtmIinkZcuCMEuI4OejDYjZRbs9LGX1BSBt/0GDAG8j3MCYNMoVcZHT1ezEMsXqE4kKs9YlFhL3IMLSmJ8OIAkGYSAJirU84IuxFSN+Qn6BY7UKR0C3W+oQjwl6EGFrLRKpQFAQNLdZ6HhBhL1L6PAECkySeXyheBrwBsdbzgAh7kaLF1y4UAW6fWOv5QIS9iOn3BCRCRihYIuUxhIlHhL2I0VozKBaRUKDIdzN/iLAXOTIxJRQqg/LdzBsi7EXOkC8ok6hCwRE0NB6/fC/zhQh7CTDoFT+mUFiEOiLJ/E++EGEvAfq9Eh0jFBZuMTbyigh7CeALGPgC8torFAaGoRnyi7DnExH2EkEmUYVCYcgfFDdMnhFhLxEGPCLsQmEg1nr+EWEvEQKGgTcgN5SQfzwi7HlHhL2E8PjEz15KFOO8SdDQRTnuUkNa8ZQQbn+AKqz5HoaQAR5/kD6PnyFfkKChOa3aOSHNj7OFWOuFgVjsJYTHb0jtmCJGa01nv5cBT2C43n5nv7eoJiJF2AsDEfYSQmsJMytm+jwB/COyiP1Bg5ODxVN7X75/hYEIe4nhlmp6RYlhJG+e0jvkLwpL2BD/esEgwl5iSJnU4qRnjHaHJ4rAavdIVFbBIMJeYgQMyUItNvxBg94xmqZ4/cE4N02hIUZF4SDCXoLIDVZc9A6l1uy50N1sHjEoCgYR9hLE7Zcs1GJBa51y3fJCbjNnGBpvEcwDTBZE2EsQj98oqhC5yYw7HK+eCoUczir+9cJChL0E0VqaHBQL4+kypLXGXaBWsVe+bwWFCHuJIvHEhY9haAZH8Zs/sLGN1rbumGV/2n2ce/64N9dDGzde8a8XFCLsJUoxxD1PdgbG6DI0v76CNU/uHBb31rZu7nx0G4sbqiZqiCkjBegKC6kVU6J4AyF/rMmk8j0UIQljlVpuaa7hG1cvYM2TO1m5ZAaPbznCN65eQEtzzQSNMDUCQSPleQJhYhCLvUTRWsvr8Qi01hzt9XC4Z4gety+vVqY/aKT0VtXSXMPKJTP41fo2Vi6ZQUtzDYMFFh0j37PCQ4S9hBE/eyxdAz7cvgAef5CTgz4O93jyJu6pTpq2tnXz+JYj3Ly8mce3HKG1rbvg8hRE2AsPccWUMCLsp+h1++n3xGZ3aq3pGvDRUO2c8PGMNmkaobWtmzVP7hx2v5zTVD38/xlVDizmwrDLxL9eeBTGN0PICb5A4cY9TyRuX4ATg96E67z+YNLiW7kiEDRSSubZfbQ/xqce8bnvPtpfUFayhDoWHlmx2JVS/wNcDRzXWi/KxjGFzNFa4wkEcdkm94vZWHVYut1+XDYLNsvE2DmpWOsANy1rjlvW0lxDS3MN3oBBmT3bIxs/voCBIclwBUe2vsm/BK7M0rGELDLZE5UCQWNMn3TIJZPYos8F2SgNUCjhrOKGKUyyIuxa6z8BJ7NxLCG7THY/+0CKk5Qef3BCxDJo6KxMfnoDhVE2opBcQsIpJszHrpT6tFJqk1JqU2dn50SddtLj9QcntZ+9f4xY8Wj6xnDZZINsFfIqlHDWQhiDEM+ECbvW+uda66Va66V1dXUTdVqByVugyTPOGuaDviCBHNc8H/Rm77PIt6hqLR2TChWJipkEFFrc80QxHmsdQkI13n3Gg2Fktydtvv3bheIOEuIRYZ8ETEY/u9Y6LbdHv2f0+i2ZMOQPJj12ooJfrW3dPLCxLel2kTDDV/d25aUwWL7fGITkZEXYlVK/Bl4D5iul2pVSf5GN4wrZwReYfLU8BsdR5zyagGGkPOE6XkbLNo0u+GXydHPilV/ieOKzLLfvT7qdP2jw8ttd3LG2lcWNE18YTNwwhUtWApy11h/LxnGE3OHxBymzT5549kwmKfs8ASoc1iyOJvIGEf/m9MDGNubXV9DSXMO3rpxJ5RO30chmZhFEowhsOUD7gufRtnIgvjDYk1uP8JOPn8sFc2uzOt5UyLcrSEiOuGImCYUS9zxRZJIN6fUHsy5abl8wYSJPtAV+ycmHuYzXuTfwfr592o85fP0jWPrbmfrKmph9oguDffi8xryIutYaf3ByvQUWEyLsk4TJ5Gf3B41xRcMkYqySuuMlWUXGiAV+9xMbcWz8D9YFz2PjvL/jwSPTWR+YR2/Laip33I+z7aXhfaILgz20qZ1X93ZldaypIBOnhY0I+yRhMvnZs/F2MuDN3iSq1hr3KGGOLc013FX1NA5jiJebP8udH1ww7G55ccan8NWcQd0fvojJ2xtTGOy2C2fzD9cs5I61rRMu7r4ch4UKmSHCPomYLO6YbLydBMdoWzcehvyJ3TARdu3ewXt6fsf6iiv43eEqWtu6hy35nZ0+jl/+I8zuTqo2/zSuMNiSpmp+9NFz2Nrem5WxpooU/ipsRNgnEZPFHTNSdFINJRxJttwxo0XZtLZ143n+u1jM0Hj9Pw5b6hFxv2lZM75pSxhqvIjyvU9y07ua4jootTRXs3rF3KyMNVVk4rSwEWGfREyGRKVAAv96ot6ha57cyfz6ilGP5fYFMs5EHcsNc+zALq7nJQbO/iSBysaY0rzRDM65EmvvAawn98QdY6LDDmXitPARYZ9E+CdBb8pEbyXRIYL3vrI/pnnFWGQa0z6WG+Ym1wYUBj1LPh0z3pEle92z3wdA2f6n444x0YlCMnFa+IiwTzKyVYSqUIkuUxztgokOEZxbV5ZyQ+hMSwyM9WAof/sJhmYsI1jRMOp2wbJ6PNPPpWzfs3HrJtpil4nTwkeEfZJR6n726AniaBdMa1s3j2zuwG4xsftof5zPPRn+oJH2wzBo6FGLfllP7MZ2cheD865N6XiDcz6AvXMr5v6OmOXGBBfjkonTwkeEfZLh8ZXuTTnSvx5xwXzjse187dFtAHzn+kWsufasGJ/7WIzVgSkZ/R5/QpdF5E2i/J3H0crE4NyrUprMdc95PwBl+56JWzeRk5kycVr4iLBPMgKGUbI3pieB1drSXMP86RV4AwYfOrdhuLVcognKZAz50stETebGmV9fwZondmDZ+ShDDRewqcuS0mSuv3ouvpozKNsfL+wTZbHLxGlxIMI+CSnV6JhEcfqtbd3s7Rrk5uXNPPFmByde+SW2rh0JJyhHY7xW+6A3kDT7taW5hh9eDBWDB3lOXTiuydzBOVfiOLwBkye2YdlETaDKxGlxIMI+CSlVP/tIcWtt6+auR7exalkTt104mwfnPMW73ryLxt9cgen/VuLa9yzo1ARx0Du+Jhx9ntEfBOcNvEQQM//w9lxWLpkxpqhH3DeDc65E6SCu/eti3DcTZbHLxGlxIMI+CfH4S9PqGiluu4/2c9uFs1i78RADL/4bZ+z7X1qnXc8PgquY5m+n/unbmbL+eykdW2tNX4oRMr7A6A20H9hwELX9EV5hMSuXL+TxLUd4aNOhUX3skYngDUNNBMrq8e56LsZ9M1ETqJMle7nYEWGfhGid3U4+hYA3EN/E4qZlzQQNzT/O2cHiHf/Ctqr3cPvxj2Jc8HnuXvRb+uffQNWb/4n1xO6UztHv8afUP3Ysa73JvYMq31H65lzDbRfOZtWyJu754z7MJpV0n+FY/N/vYqf1LCyHX49z30zE3IlExBQHIuyTlFLzsyezVs+t6OEDe9ew27GYG47dytLZtazdeIgzZtRw4sJ/wLCWU/vnuyCFN5igoeka8I66jccfHDP2fU73nwgqM9/bN4d7X9nP2o2HWL1izpjJY5FY/Ic7G5jBCZbWuGPW59piNwydcdVMYWIQYZ+kJGr6UMwkmzw8v+85LBjc2vOXnNFQy/M7j7NqWajeiuGcysnlX8HZ8Rqd69fG7Jcs/HDAG0g6kRoIGhzv847p5nqXdwO+087n0nNO51fr21i5ZAY3Lm0aczI3Uq53+lmXAPDmq8/FrH9174mctsiTVnjFgwj7JMUfNDKug1JIJLRWtca24yHWs4hFCxbwVkcfly+YxtqNh4Zj2PsXfpzuqoXM3vxd3tp7CBi7lszJQV+cr1lrzfF+LwFj9L+ppa8N28nd7Km6aLim+uNbjowZUx9drvf9l12Oz+TgxK6XeWjTqTH/v99ty2mLvFINky1FRNgnMe4S8rNHhD26jID9yOuUudt5p/5q/vx2Fzcvb2bjgW5WLWs6FcNuMuO+4vvU0cPxZ/45pVoyWms6+70M+YJ4wt2WTiQQ+0S49q8D4Gs7moZrqkdXdExGTLles5XgjBaurGrj3lcOxIz5vJmplUpIB4/414sGEfZJTKn42X0BY7jQVnQZgYrdv8VncvK9A/O47cJZwyK6duOhGGvcO70F9+z3s8r8Bx5c/05K4Yf+oMGR3iEO9wzR0T1EX4px7q4Dz9HlmMXt11w6fI5UEqZuWtYcMyZP/VJmDO1h1bm1w+6cluaanLpLxGIvHkTYJzFDvvhIkmIkWnAiIvm9J7Zg3fUYTwXexS0rFnLj0qaY9SNFdPP0GygL9vKdM95OyTWSDsrbh/PweiwLPpigpvr4EqY89UtRRoD9W1+JcefkKmplMlQGLSVE2CcxhtYl8Xo90r/e0lzDF2ftxWkM0D3v+mFRj14fLaKtbd3csb6K/rLZXOV5KiXXSDq42l5CGQEGZ12R8bE2BUKNNb545skYd85rOWqRJxOnxYUIe7ES9FO++2GmvPpt6p7/PPWPr6Lmte9i8vSM6zClUMY3UcbpjAOP0Wut49/3zhhToHcf7ecb1yzE23IbjuOtnO9oG1ctmVRxHVhH0FGDt/68jI/11kkzfeWzme3ZDpx6E9nW0ZtSrP14kcSk4kKEvdjQBmV7HqVp7XuY9vznqdryCxwd6zEPnaB6809ouu9Cqjb/FBUYSulwpRD2GG2xt7Z18x9PrOditQV99kf4f9csGtP6jviv+8+8EcPiouqt/x23a2RMjACug3/APfMyMJkzPtxNy5qh6XwcRzYNx+BHxpwL61os9uJChL2IsJ7YTcODVzJ93R0YVidHr/ol+1fv49AnN9Lx0Wfp+OizeKefy9TXvk3946sgOPaEXrGHPUZPnELI+v5+SxcmHWDg9GvGVclR2ysZmP8hyt5+LK7IVqY4jm7C7O3BPet9KW1vt5qpdtlG3cZbvxSztwdrT2zserYnOfUE13sXMkeEvRjQmort99Pw2w9iHjzGsSt+QsdHn8M96wpQp9LQfbVncfSaX3H8vT/CeWQjUzb8c8LDjWzu7PYHeXVvV06TW3LFyKJUNy1r5kzvWxi2CnxTFwLjm5jsPftWTEEPFTsfzOo4y/Y+hWG2425+z+jb2S00TXHRUO1kSpltVHH31C8FwH70jZjl2baupaJj8SHCXuAo3wDT1v0VdS99GU/9u+j46DoGz7gOVPKPbuDMG+k762aqW3+G68C6uPUjmzv/eU8nd6xtzWlyS65IZEk6jmzEM+Ndabk8/FMXMDRjGRXb70upzEBKaIOyfU8x1LQCbStPuplSiqllNqzmU59tjcuK3Zr4Ovw1cwnaq3Ec3RSzPNuRMZk8KMz9HTg6XsPR/jLOQ3/C0nsgewMTkmLJ9wCE5FhP7mH605/C2rufk+d/hZ7z7hhV0KM5cdE3sR9rpe75L9DxkWcJVDYOr4tu7rxyyQwe33KEn338XC6YW5urS8kZI90OJncXtu636Z//4bSP2b9wFdNe+AKOw+vxNLw70yFiP/YmloEjnFz+1VG3q3BYsJhjP1+lFHXldjp6huKtZmXCO70F+7E3YxYHjFAnKas5O3ZbOhOn5oHD1Lz+Iyp2PoDSsft7py5gcO5VDMy7jkD17KyMUYhFLPYCpeztx2h46CpM3l6OXPsbepZ+PmVRB9AWB8fefw/KCFD74hfj1kc3d04lIadQGWmxO45sBMBz2vlpH3Nw7tUEbZVU7Fg79sYpULb392iTNeQ6S4JSimqnNeE6m8XElCQuGe+0xdhO7o6bLM+WO0ZrPb5ENq2p2fAvNN13ERW7HqRv0Sc5vPIBDl/3Ww5f/whdF30Tw1ZBzcZ/pWntJUx77q+wntyTlbEKpxBhLzCUb4C6F/6O6c99Dl/tQjo+8gyehgvSOlagejbdy/4OV/vL2I9ujlkXKSgVSW7509ud2Rj+hBI0dFzSjPPIBgyzA++0JWkfV1udDJxxPWV7fz/u8NH4g2nK9v6eocaLMezJXV2VCaz1mPVOCyYVX9bXW7cEpYPYunbELs9SeKJ3xOT0WFRt/gk1m/4N95wrOfTxP3Pikn/E03QxnoZ34zntfPqW/CVHPvQobbduovec1bj2P0fjry9j2rOrsXa/k5UxCyLsBYX9yOs0/uYKync/RPd5f83h6x4iWD4jo2P2Lfw4QXsVVW/+bHhZdEGpSHLLV367lVdzlNySKxL61w9vwFvfAubRI0rGon/hKkxBL+V7HsnoOLaubVj7DzE496qk25iUGjMCRimFyxbva/dOWwyA/fiW2OVZstjHY6279j7F1PXfZWDedRy/4icEKpuSbhssq+fkBXfRdssGes67A9eBF2j89WXU/uFLmPs7sjH0SY0IewFgHjxO7Ytf5rRHPwTA4esfoXv5VzMWJwBtK6dv0S2U7X0aS88+YERBKUJuma9fvYA3D2VonU4wI4Vd+fqxdW3HM2N55seuW4S3bjHGpl/SejA29DFZSd9ElL3zJFqZGZyTPMyxwmEZtclGBJc9fkosWFZPwDUtobBnI5Il1YYstuNbmfb8X+OZfi6dl/1rTLTWaBjOKXQv/yptN79G3+LbqNj9MM33XUjti3+PpWd/JkOf1Iiw5xHl7aN6049ouj/sj1x8O+0fXYd3xruyep6+xX+BNtuofvM/gfiCUhAS95uXz8zqeXONNxgrOo4jm1DawHPasqwcv2/hKqYN7eWRJ58YjiAaq6RvDBE3TMMFGI4pSTcrSyDYiXBZzaiRgqkU3mmLsR/fOuLUOmOr3TBSO4byDVL/9O0EHVM59sH/QVsc4z+Xq5YTF32Ltk+8TN/CVZTvfjjkg3/2szjaX0m5N60QQoQ9D9g636L2xb9n5i/PZcqGHzDUtIJDH3uRExd9a9RwuHQJuuoYOPNGync9hNmd3JdebO3y7n1lf0w8vuPwegxl4b726Vk5/sAZ12FYnPxT8ybWPLlz1JK+I3MDAPZu34itd/+obhiLyYQjSTjjSEwmhTPBtt66JVi730b5BmOWZ5pVPORPrUhc1db/xjJwhOPv+wlBV11G5wxWNHBixXc4dPNr9C75NM62lzjtsY/QdN/FVG/6UWguQWLqxyQrwq6UulIptVsp9Y5SavSYrsmG1pjdnbj2P8fUP95J068uoPHBKynf8wgD866j/canOfaB/8p52FfPOZ9BBX1Ubv2fpNt4fMVjFWmtmTctNh4/uP8VtunZzGmYlp1z2CoYmH8DMzt+z6qFtrgIomgxj+QGRJpSt7Z1c+SP/4OhLLjnXJn0HGX28cXaJ9reO20xCo29a1vM8kzrAKXyoDd5uqlq/RmDs96X1TfNYNl0Tl74ddpu3czxy/+dQMUMpmz4AY2/uYLm/11K7R++RMX2+7F1vgVBX9bOWypkHMeulDIDPwGuANqB15VSj2utd4y+Zxr4PSjfACqnKfBha0BrFDpsHYR/DCMUk6uDqKAfFfShgl5UYAiTbwCTrx+TtxvL4DHMg8ew9rVhO7ELczg93bC4GGq8kJ6WzzB4+rUYjuocXkcsgeo5uOdcSeX2++he9kUwxX/0AcPAGwhit2ReyyTX+IIG5zRVD8fjf+jsKdzZ/RYn534iq6GbPed8hort91O7/V5uXv4lHt9yhHOaqmlprhkW84gFH2lKffmCaTzx+tv8yfIS7tlXj2rFpuqGieCyWVDKF2NJ++pCE6idu9dTERXmuWHfCR7vO8znLj19nFcdIpWJ06rWezD5+uk+/+/TOsdYaKuTgfk3MDD/BsyDR3EdfAln2wuU7XuKyp2/Dm1jshIob8Bf2USgspmgq5agYwqGowbDVoFhdaEtLrTFjjbZ0GYrKAvaZAZlRitTKJRYKSDs6gq7vDRRrq8YN1hqcwiH0BnKAAAgAElEQVRxWGzgqARzblOIsnH0ZcA7Wut9AEqpB4BrgewL+3N30fj6L7J+2GyjlYmgaxqB8tMYnHMlviln4qs7C8/0c7MyIZou/fM/TNm+p3G2v8JQ84qE2wz5ikPY/cFTha9WLpnB9g3PYrMHqJy/AvcY+46H1/tq2Mlybra8QMe7vsc5TdUxYj4y0evyBdNYt/M4P5m7CWvHAMeXfCrpscfjholgNikcVlOM6AbLpuF2TKd9x6t45t5CS3PN8FzADz+SXtinP2iM2bjaPHicqq2/YGDedfhqF6Z1nvEQLKunf+FN9C+8CbTG0ncQ+/Et2Lp2YO1rw9J3CNf+ZzF7TqIK2Sf/8Ydh3uU5PUU2hL0BOBT1/3YgLjtEKfVp4NMAzc1pVs1bcA099gaCY/SVzJjhJ7NCDz/FFZgiT3cz2mxFm+1osw1tdmDYKzBslRj2KoLO2qxU8Ms2Q83vwbBVUP72Y8mF3R9k4t4j0icSEROJx//e7JNwGDYF5jCaxNitZvzjiM3efbSf8y/+ArY/f5TKbffRcu5nh4uKtTTXxCR6XbFgGhsPdHPL+Q0s3vIbumrOwTu9JemxXeN0w0Qos1virGk9o4VLju/iovBD5qFN7dx24SzObjgVO//q3i62tveyesXcMc+Rihum+o27UUYg9AY40ShFoGoWgapZDM67NnadNjB5e0MC7xvE5B9E+QdDb9iGP/SWbRigg6ADKB1+Mx/OkA19N1TMdyTq9wx8/E6bBVdtem9Q42HCSgporX8O/Bxg6dKl6f1l5ryH/uplY1oSQiwPbGxjfn0FLc01DM5+P679z7Bl/1fY2emNK47l8YfC5OKiLwoMX8CIice/aNda3I7p3LXuON+wT03ojnHazNRXOggampNuHwOesX3Qob9PM+79F1G15b/oXXL7sKDDqQfLFQum8fzO46xeMYdPTtlO/ZZjfLH741zQ1p3UNVRmS+/2cyWaQJ22mCn7n+HGsyv4Rfghs3bjIc6YXsH15zby6t4u7ljbyo9XJX/QROP2ji7s5sFjVG6/n/4zP1p4ZQGUKeSGcRReNrVy2XCV5f6tPRuTpx1AdCZCY3iZUCBEF/0amHctZm8vLz71m4Qhe1rrooiO8QWMmHh8e9db6PolSUv0OqwhUVdKYTGbmFbhYEaVM+nxR0a59J77V1jcx9j5zH8NL4t+sMyuLWP1ijms3XgI84af4a9o5MKrbklaLthsUjgTJBylgsVs4sFNh2LGF8m0ffvNV2Kadn/ziR18/5ldw6KeSj0gf9AYc+K1YudvUIafnpbVaV2DkFuyIeyvA/OUUrOVUjbgJuDxLBxXyBLRvuB7DjXRQzlfatiW1JIs9CbXQUMTMIzheHzlG8TavRdv3dkJS/TaLKZhUY/GaTNT7khsNY+sgPmqsYjtzOGyo/+NyR3K0I1+sNy0rJkblzbxw4s0dSc30bf4ds6ZVZe0XLArTWs9wuLGqpjxPdQxFYBbZ3XHNO1eNquGn760l0+c35xykbcx32S0QcXOXzPU8G4C1XMyug4hN2Qs7FrrAHAH8CywE3hQa7090+MK2SXiC/7fDYd5e+qlNB5/KWmXpUK32EdmnNpO7ECh8dUtSrh9jcuGKUlm5xSXLaHbKfpheO8r+1nz+10cueT72P19TH/us2AE4hK9TN5eLt7+dYL2KvoWfGzUa0hUHmA8XDSvLmZ8P93YQ499Bu+yHRge/6plTfz57S5uvWAm921oS7lkRP8Ywu5ofwVrXxt9Cz+e0TUIuSMrcexa66e01mdoredqrb+djWMK2SW66NfPT7Rg8g/iOviHhNv6AoXdkX6ksNs73wLAW3d23LamJDVWIljMJqqSVFUcWQFzztnvpus938PZ8SpT1o9oYhL0hkos9+zj2JU/R9srk55TqcSJRuPBYTHHjO/GpY1YGs8bLi3Q2tbN2o2H+Pb1i7jl3bO4+6ZzuGNt65ji7vYFCIwRnFC5Yy1BezXuOR/I6BqE3CGZp5OAkUW/rrz6Brqowr/l4aT7FLLVPrKUgL1zG0HnVIJl9XHbltktY04EVzutCWu1jKyA2drWHdXE5KdUbPs/LL0HIeil7g9fwtnxKp2X/SuexotGPZ/Dakr6BpEqNouJLe09MePbZz8Ta/8hzO7OuHpAC2dU8uNVLWxt7x31uGNZ66ahk5Tte4aB+TekVTpAmBik0cYkIK7o18xaemd9gOa2R2kPeBLeoG5fgPJxJs9MFJEY9gi2zrdC1noCAa9I4kOPxmQKVVc8MeAdXhb9MGxpromNX7/4W9i6tlH3xztjjnNy+VcZmH/DmOdzWTP/u766t4tvPb4jZnw/fGInvwTsx1q5aVls0bE+T4Bls6aM6mcPBI0xyxBU7H4IZfjoW7gq42sQckdh3rlCVkk0gedadBXWAw/gOLyeoQR9OAu1vMDIxsoq4MHWvYeemZfFbWs1p54AVOmw0OP2DbugElXAjI5fP3LdQ9iPbMLa34Glv51AeT39Kfqc042GiWZrey8/+PBi5k4rHx7v2UsvIfiGGfuxzbhnv4/Wtm52H+3npmXNaK3pGfJTW25PeswBb2D02jBaU7FjLZ7p5+KfembG1yDkDhH2SYqn4d0YZgeug39IKOwBw8AXMLBZCstb5wvGlqO1ntyNMgL4EvjXx5Our5Si3G6hd8gPJH4YRseva4sz1EBinOO3mk1Z+ZuuXjEXjz/I4Z7QBHgoiqed2yvn4zi6OeaNI0K/J0C105qwoYc3EKTH7R/1nPZjrdi636Hz0n/JePxCbimsu1aYMLTFiafxApxJJlChMP3s45k4Ha8rqTLJJGo2yYa1HsFuMQ3PH0TeKJ7rbUQdaeWfntgWV4VSa013AvEOGprjfd4xs3HL33kMbbIxMEq1SqEwEGGfxLibL8PWu3+4AcdI0mlinGsSCXvQVkmgMj52fbyWsdVsyqrwJiLTMMdolFLYo66xpbkGS/O7cBhubp/vS5inMOANxHyuWmuO9XnGzuY2gpS9/QTumZeNGvEjFAYi7JMYd9gv7Tr4YsL1hZio5BshQLbObaH49RETp+lO/FY6cme1ZyPMcSQPvH4opgnI/7aFKkke3fFyXH14CAn54Z4h2k646ez30tnvTekB7jiyEYv7GAPzVmZ1/EJuEGGfxASqZuKrnourLbE7xtC64Kz2GIs96Md2YmdCN8x4qyZGKLNbsJhyc1s4E3VAypCWpqrhOvBrntzJhcuW0UcZN9Yfi8lMHUnAMOj3+BnwxoY3JmoY0trWzdFX7sewOHHPuiKr4xdygwj7JMc98zIcHa+h/ImzUAtJ2APB2MQpW/fbmILeOGEf6aIYL6mESKZDLtw8l5wxjW9cvYB7XznAslk1rH29A3fdOSww9iStmzMaI0sptLZ18+0ntnFW70u4Z78PbXVl/RqE7CPCPskZmnkZpqAXR8crCddn2l4tm8S5YcIdg3y1saUEoicV06HCMXZSUzqU5UDYzSbFstlTuXFpI+t2HmflkhnYZy7DdnI359bbktaqSUZcKYUnd3L3+b3Yfd0MnH7t2AcQCgIR9knO0GnnY1icScsLePxBjAIpLxBXI6ZrB4bZgX9EIap03TARLGZTVic5IVQHPlGYYTZ4a0QG6g7TPJQ2sHduSet4Lc01zK0tGy6lcE7vHwjaKnlNncMDG9uyPHohF4iwT3bMdoYaLw4Je5JwN3eBuGPiImK6duCbOj+uqYnDmvnXOtvumPIMqzkm49W9Xdz1u23D5SK+cfUC7twYyiS2H2tN65itbd3sPtaP3WLiic0HsL/zFIemX8Y3n96bsNSzUHiIsAu4Z74Ha/8hrL2Jwx4LJTrGGy3sWmM7sQPf1Ph+SY4stPZz2SxYs2hhp9staSy2tvfybx89JyZD9m+uWcZJeyOOo5tjtk02MRpthUcSm9ZcexbfuX4RF/EmtkA/3z24IC4uXihcRNgFhhovAcB56M+J1xeAsBuGjom1NruPYfZ045u6IGY7myXzAlsRRrPaUxHJ6DFl8yERzeoVc7lkXmyz7JbmGqyzluM4sgGien8mmhhd8+TOGCs8upRCS3MNn5nSSpeu5MT0C0TUiwgR9klKtDAFqmbhr2jCu+eFhMIUMAy8gfyKe/zE6c7Q8hFNlDP1r0dTPkplyFREMkK6LfBSxWRScclYQ00rMHu6sYUzcyHxxOhIKzy6xvxb+9o5vefP7K17L293eZKGTgqFhwj7JCVGmJSivWYZlUdf48xpidvF5dtq9yZorgHEWezZFPbRJlFTEckI46lZky7Rwv7AxjY2mhYD4GoLJZ9F3iZG1phPZoW3tnWz4en7cOKj8ZJPDl+riHtxIMI+SRkpTD9ta6aCIc53HEy4fb7rxoycOD2yexNuRz2Go3p4WWtbN//36oGsnne0TNRURDJbRb/Gwh41rzC/voK71h2nu2ohzrY/xrxNJKoxn4jdR/v562lv4i9vwFt/Xkx1S6HwEWGfxEQLU82iy9GopH52j9/Ia9jjSFdMk38/mzynxblCzmmuTrR72jht5qTCnIpIToS1DgwnZEVcad+4egG/7Z2P7cgbfP+x17n49FBP1OiGK6NZ4avOLqe+6zUG510LKnTsRP1khcJEhH0SEy1Mv9nupqdqQVJh11rnzWofWYOdoJfKwf1MO/28GFfIdz+0KOWGzeOh2mWLWzayK1UykZyoZiWRpKyIiw1Anf5ezAR5l36LS8+cNmqN+ZGU7X0SZQQYOOO6CRm/kF2kHvskJVGHoEeeOIPb+p5E+QbQtvK4fQZ9gQmzQKPxBmJrsNtOvo0yAtTMaWFleeiN4+blzVx4et0oR0mfcruFbrMpJipnrEYcELLWJ6qevVKhCdTIOL7x2HYwqrnJ7OQSUyhRaawa89GU7/kdvpp5CcNJhcJHLPYSp9xuoaHGSaXTiikqwiORMM1bfjUmHcB5+LWEx3J7g6N32MkRXn/iidMtvsaU/MXZoNoV62uPjh6JMNJVkaxJdq6Iro8TMDSDARNt1Uu5yrWTNU/sSPnvY+7vwHlkQ8haz0FpBSH3iLCXMFazidpyO3aLmdpyOzOnurCHo0YSCVPTkksxzI6k7hhD67zUjvmvP++LESVb1078ys7fPt8f4wr5uwe38OrerpyMocJhHVcsusNqzmqETipEhP0Pu45jMSluXt7MI71nUjZ0mO9f6kx54rNyx1o0ioF51+dyuEIOEWEvUUxKMb3SEZOso5SiehQrUlsceE47H+ehPyXdZtA7ehf7XHD6tPIY/7WvYyu7jAZuuXDu8MNp2eyp/HhVC1vbe3M2jipX6hb4RFvrEAp5bG3r5uV3TrDm2rO47cLZnPfeUHPt+uMvpzTxqQIeKrf9H+5ZlxOompnrIQs5QoS9RKmtsCf075bZR0+VH2q6BFv325gHjiRc7/ZNrDsmEDRY3Fh1KjTz5X1YurZT3ryEG5c2DW9nt5q4YG4tq1fMzdlYKlL0mVvNprzMRdgtZvYcG4hxsc0/82z6y2ZR1p78YR1N+Z5HMXtO0rvkL3M5VCHHiLCXIA6redRojNF6ew41RcoLJBYCQ2sGJ9AdE0lMioRmPrNhK1Pop2JmS8x2dnPu3R4q/BZkHqNkwUh//ETyyQtmxbnYgvPez9z+TZj7O0bfWWuqtvwC79QFeBouyOEohVwjwl6CjNWUudJhSSpOvqkLCDjrcBWIOybS6CMSmnnHQg8A24ONMdtNVPSJ1WxiWoUj6fpql42KHLbXG4tEfv3exX8BQHXrz0bft/1lbCd3hax1mTQtakTYSwyzSY3Z0EEplVx8lGKo6WKc7X+OKSAVjds3cTXavQEjJjTz6uknAfj6a8RMqGbSMWm8OG1mppbZ45ZXu2xMKYuPeZ9IEv0dghUN9J95IxU71mIePJ5036qtvyDgrA0lJQlFjQh7iVHhsKbU/adylC5BQ80rMA+dwNa1I+F6rTWDvtxb7VprvAEjJjTT3rUNf3kDf3PNsuEoD6s5exUdU6XKZWVGlZMpZTbK7BZqCkDUIfkDrufcv0IZfqre/M+E663d71B24Hn6F92CtiR/IxGKAxH2EiPVBhGjFbgaarwYAOehPybdv9+Te2GPJCZFh2baurbjqz0rJmZ8otwwI3HazFS7bEyvdFBTAKIOoc81UTPuQNUsBuZdR+W2/8M0dDJmnfL2Mf2ZzxC0VdK36JaJGqqQQ0TYS4jxNodI1lw5WDYd79QFo4Y9evzBnDe6HlnRUfndWLv34quL73EqnOLBTYf41+d2x7iqWtu6udtzNSowRNWW/zq1cdDP9Gc/g7XnHY5d+XOCrtxk7woTi5QUKCEqneP7OF2jJNAMNV1C1dZ7Uf4htDVxKd8et5/6qtxFo4ysAW87sQOFxjuieXW+LPZCZUljNfe+sp+Xdney5tqzAEIlBijj1oYraNj8Y+yd2+g/88M421/GdehPHL/sh3iaLs7vwIWsIXdEiWAxmXCNs6GDZZSSskNNK1CGD8fh9Un3d/sCOW3AEVdKoGs7AL7as2KX56g7UbFy8Rm1w4L+tUe38bVHtwGw5tqzCHzwX+lp+Ry2EzuY/tznqNyxlu7zPs/Ago/mc8hClpE7okQoS7OnZrKHgee0ZRhm+6h+doBetz+t846FL2DEFN0CsHduJ2ivJlDRMLzMYjJhEWGPwW4xce7MKXzo3Aa8AQNvwOBD5zbQ0lyD4ZhC97vvpO2WjRxZ+Ws63/PPdJ//9+Nq9ScUPhndEUqpG5VS25VShlJqabYGJYyfdDMdk02gaosTz4xluJLUjYkw4A3ECXA2cCeIurF1bQtZ61HRPOKGiUcpxVvtPTyyuQO7xYTdYuKRzR2xwm0yM9R0Cf1nfQKUaVyt/oTCJ9O7YhvwISC1fGUhJ1jNprQLTtktppiqj9EMNa3AdnJX0vICEXpyYLXHZbcaAWwnduGtE//6WLy6t4s7Hw31Ov3O9Yv4zvWhv9k3HtuetMLjeFr9CYVPRneF1nqn1np3tgYjpEcmdUmUUkmtdvfMS4FTfTOT0e/xZzVCJhA08I44nrX7HUxBb7x/XYQ9jq3tvXxgUT1rrj1ruN76mmvP4tL5daNWeEy1H6pQ+ExYVIxS6tPApwGam6W9VjZJJsyp4rSZGUhQJsA/ZT7+ikZcB56nf+GqUY/R2e+lscaZUnLUWLgTPCTsXaEJwLiIGPGvx7F6xVx8AYP2bvfwsmQNNaIZ2ervnKZqEfciZcy7Qin1vFJqW4KfceUda61/rrVeqrVeWlcnsbLZIhM3TISk0TRK4Z51Oc5Df0IFPKMewx806M6SS8btjRd2W+c2DLMDf82p6o2mcNcgIR6bxZS0HtBXH97KQ5sOxSz70bo9fOXht8Zs9VfptDK90sGMKiczqpwJk6GE/DPmp6K1vlxrvSjBz2MTMUBhdLJRHtZsUsMNOEbinvleTIEhHB2JuypF0zvkzzj80TAS91a1d23HN3UBmE5dr1VEfVSSPfDPm1nDPX/cNyzuD206xONbj/DBRfWj9kMtt1uoLbdTZrfgtJlx2szMqHaMKylOmBgkQanISTfMcSROqznOrw3gabgAw+LEdfB5hsI+92Rorens99JQnb5LZsifoN671ti6tjM49+qYxeKGGR2nzZywEmekjv09f9zHy+90sa2jj8+umBNT3x5i3Td2q5m6ivjCZ1aziRlVDo70enISHSWkR6bhjtcrpdqBdwO/V0o9m51hCalgNZuwW7In7InQFgdDjRfjOvA8pNBgwxcw6Oz3pj2ORMXFLP3tmL29EhEzTpJ9phAS90UNlbzV0ceihso4UY/GYjIxvcKe9GFtMZuor3Ikja4SJp5Mo2Ie1Vo3aq3tWuvpWuv3Z2tgwthks0uPw2pKeuO6Z12Otb8d68k9KR1rwBugx+0b9xgMQzOUoInHqYxTqREzHqxmU1I3yUObDrGto4+zGyrZ1tEX53OPZmq5bcwkMKvZlNcGI0IscmcUMdlyw0Ao7NFhjf86PLCxjY2W8wBwHXweSC0j8eSgL2GS0Wj0DPkJJqjzbu98C61MIR97FOKKGZtEhd4e2nSIe/64j9Ur5nD3TS2sXjEnxucejc2Sepu/Kuf4Gn4LuUM+hSIlm26YCIle3efXV/C1F07SXbkA14Hnx5WReLzPm3J8uz9o0DuUOKrGfmwzvilnxhQjy0cN9mIk0Wf6xsFuVkf51G9c2sTqFXN442B88lKNK/VyxEopppYXRvniyY4Ie5GSi2bJiaIoItERD/YtxH5kE3c/sTHljERDa472elIS95ODPn694WB8vZKDXZgOv4G3/ryY5eJfTw2n1RznYvveDYvjfOo3Lm3iezcsjlk2Hms9gstmyUsjbyEWuTuKlGy6YSI4rOaEE2AtzTUYZ3wAEwZ/27RnXEkrqYj7kC/IoDeQsF7J2ifX4QgO4hkp7PLKnxImk0p7LqJ6HNZ6NFPKbFlJVBPSR+6OIiQXbpgIiXyyrW3d3LOnkpO205h+6Omk9UaSERH3Xrc/LpQxaGi6BkJRNInqldy1OBRH7amPrTEnFnvqjBYdkwyr2UR5mpZ3JvsK2UHujiIklzfNSHfMcCPpaxaiFn2I5Wob//HE+rTE/cSgl/buIQa8AQa9AY71eWg76Y6Jfx5Zr2SebydBxxQCVbNijifCnjrJOmWNRqat/qpdqfXeFXKD3B1FiCsHbpgII6276EbSg/NWYtJBfnj2wVGLSY2GP2hwvM/DsT4Pg94AWuuYWuCReiVXLJjGQ5va4dDGkLUeJRImpST6Yhw4rOak5QUSkQ2LW6z2/CJ3R5GRSzcMhCzh6PofkUbSD2xsY8PgDHzVc1lw8nluWtactUYMEd/6Q5sOsebJnaxa1sTGA918blk1FYMH2OdYGLO9lBIYP+MR2WzFo9eI1Z435A4pMiocubeCHLb4r8X8+grW/H4X70y7AkfHenbu2cOaJ3dyuGco4847Ed/6va8cYNmsGtZuPMQ3rl7AR2YcBWCLnhezvUycjp8qZ2oiazWbqHBkR9gtYrXnDblDioyJuFESTbZFxPeuPWeg0Ly17ld84+oFXHrmtKx03mlpruHGpY2s23l8uBa4/egbaGXmwhWxCc32BIlUwuhYzKaUIqmynT0qVnt+kMdpEeG0mSekv2eojG98vZeW5hrePGcZOzc38bGK1wmGwx4jkSwrl8zg8S1H0uq8k6gW+AePvoGv9qyYxCQQiz1dqpxWBjzJs4Fz4RePWO39ntz0xhUSI3dIEZGtV+SxMJsS1zmPiG/HaVcyc/Atdu3aAWTeeWc48iaqFvi3n9iG9WhrXPw6SI2YdLFbzMlr7xOKhMmFdS1W+8Qjd0iRYFKKsgw7JY2HkQIQLb5nXH4bAPtf+AWtbd1x1vZ4QyGjI28g7JZp6sUSHIqJX29t6+bBTYdEJDIgmatlSpktZ24+8bVPPCLsRUKZ3TKhgjay3V60+AaqZuJuuoTbHC/x0s4jcdZ2os47oxGJvInmyqrQ5Otm43Tg1INlcWNVhlc2uXFYzVQ4Yi3oSqc17SzTVBGrfWKRx2iRMBHRMNHYLSZMSmGEM0VvWhbbp7Zv0S3UP/0pLmEz77n6qoSddzLplznPs5UhxzS++kIvK0/sH/bdXzhX2ipmSl2FnRqXlT5PgKChqS2Pb6CRbcTXPrGIsBcB2ehrOl6UUkk78AC4Z11BoKyeq3zPcLT5EzHrUmmcPCpBH662lxiYezUr7afxq/Vt3Lw8ZNVLxml2sJhNTMkwu3S81LisDIST0oTcIndJEVDpzE8Dg+hU9OjsUABMFvY0XI+z7SUsvQeze97D6zH5+tlWfkGc714mTosXi9lE5QS/eU5W5C4pcMwmlbebITqePVHlxb/ZswStTFTuuD+r53UdeJ6gyc7fvl4T57tfv/9EVs8lTCw1Ltu4yhsI6SHCXuBUOvI36RTdWi1R5cXV11zE0Oz3UbHj1xBMv89pDFrjOrCOfRXn8eVrzonx3X/3+kVsbe/NznmEvGAyqZxP1Aoi7AWNSam8uWEiREfHJIpX71t0C2bPScrffiwr57Oe3I21r43qlmvj/PQXnF7L6hVzs3IeIX9UOixSxC3HyF+3gKlwWPL+2hodz54oXn2o8SK8tWdRs+luMMbX4zQRZQfWAaEG2iORidPSQFro5R65UwoUpRRVebbWITSBajaphNmha57cSeuhXrqXfQlr7wHKdz+c8flcB9bhrVtMsKw+bl0uq1oKE4u00MstIuwFSpl9YurCpILTZk6YHRqJV3fPugJv3WJqNv0IgunHKZvcXdiPbmZw9vvi1llMpry/vQjZpbbcLp9pjigM5RBiMJsUU8tynzSSKmU2S8Ls0JbmmlDiklKcXPYlrH1tVOx6MO3zuA6+gELjnnVF3Dpxw5QeZpOakOSoyYjcLQVITVlhhYS5bPGd7kcyNPMyPNNbqN50d9oRMuXvPE6gfAa+2rPi1kn8emlSZrdQLrHtWUfulgLDYTVTOUFVHFNFKTV2Q2Sl6F7291gHOqh665fjPoft+BZcbS/Rd9YtMW3whteLsJcstWV2iZLJMvLXLCCUKtxX01T6rA41XcLgzPdSs+H7WHv2juv4NZvuJmivonfxbQnXi8VeuphMivoqR0G9pRY7crcUEDUua8FapmWj1PEeRim6Lv0B2uyg7oW/BSOY0rFtXdsp2/8svYs/hbbFd16ymk0FM5Es5Aar2UR9lQOTVIDMCnK3FAgVjtyXTs0Es0mlVIgsWDadrkv+CcfRN6h68z9TOnb1prsxbBX0Lb494XpphTc5sFvMTK90SHnfLCB3TAFQZrdQV1GYLphoUrLagcF51zE454NM2fADbF07Rt3WemI3ZXufovfs2zAc1Qm3mejKlkL+cNrMzKhyYDGJNGWC/PXyjNNmZloRiDqQUjNkAJSic8V3CTqqmfHYR7EfezPhZr9dvwfb83gsPboAAAkUSURBVHehLU56z/lLIJTd+sDGtpjtxL8+uXBYzZxW7ZAHegbIHZNHqpxW6ovo1dMyjrrwhquWw9c/jGErZ8bvPozz4Isx682Dx/nU3s9T17WeNxd8CcMxZTi7dX79KT+7SSnJOJ2EWMwmZlQ5qHblpg9rqSMBpHnAYjJRV2GPqXdeLJTZLXj8qU2KBqrncPhDj1H/5Ceof+pWBs64AV/N6QRddUzZ8H1Mnm42LP0Rq984jZXmU12SohOhxGqbvCilmFJmo9Jhodvtl+5L40CEfQKxmk1UOq1U2C2YijS0q8xmZjwV0YNl0zh8/cPUvfhlnAdfoGLXbwAIlNVz+EO/Y3rdIlYG9g93Sdp9tB9gWNztFhOv7u1ia3uvVHacpFjMIUOo2mVl0Bug3xPAHzTyPayCJiNhV0r9ALgG8AF7gdu01j3ZGFgpoJTCZjHhtJpDP0VooY8k4o5J1WoH0LYKjr//ZwAoXz/W3oP4K5vR9sq4ipGrljUNFxtraa5hy6Ee/u6hLfx4VUuuLkkoEqxmE9UuG9UuG95AEI/fwOsP4g0YIvQjyNRiXwfcqbUOKKX+GbgT+Ermw0pOvjLUot18SilUeJlCYVKhZWaTwqwUFnPox2Y2laR/cDzumJFoWwW+ukUAMRUjW5prOKepmjVP7hwW95VLZvD7rUf58cdbuGBubTYvQShy7BZzaO4lXAFVa40/qAkYBkFDYxgQMAwMDRqN1oR+Ir+HjxPdf3UiWrGaJ0gPMhJ2rfVzUf9dD3w4s+GMTX2VI9enEMag3G7hxEDmHZNGqxgZaejx+ctOF1EXxiT0dqywSTwIkN2omNuBp5OtVEp9Wim1SSm1qbOzM4unFSYas0llxa2UrGLk/PoKHt9yhE9dNJv7NrTx6t6ujM8lCJOJMYVdKfW8Umpbgp9ro7a5CwgASbsaa61/rrVeqrVeWldXl53RC3kjV00Sot0zX7jiDH68qoU71raKuAvCOBjz7tRax/coi0IpdStwNfBerSfCSyUUAuU2CyeUj2x/5NHuGYfFxAVza/nxqha2tveKS0YQUiQjV4xS6krgy8BKrbU7O0MSigGTSVGWoTvmgY1ttLZ1xyybX1/B7qP92CynCn9dMFeaWAvCeMjUx/5joAJYp5R6Uyl1TxbGJBQJmbpj5tdXhPqmhsU9OvPUlWJdGkEQ4sk0Kub0bA1EKD5c4UbXQSM9d0wkCiYS2hideeoqgZh/QcgXEhskpI1SKmOrvaW5Zji0ceWSGbQ016RcIlgQhMSIsAsZUZ6hsI/MPG1t6y6JDF1ByCci7EJGOKzmtLOBo0Mbb7tw9rBbZltHb5ZHKQiTCxF2IWPStdoTZZ7+wzULhwuBCYKQHhJ6IGRMhcNCt9s37v1uWtYct+zdc6cyo8qZjWEJwqRFLHYhYyxmU9YyUV1WsTUEIVNE2IWsUOmwZuU4MnEqCJkjwi5kBact/UnUCA6rGZv0NxWEjJG7SMgalc7MrPZM9xcEIYQIu5A1KuwWTGk2ErCYTBnXnhEEIYQIu5A1TKb0M1ErHJaS7DYlCPlAhF3IKlVO67gFWiklbhhByCIi7EJWsVlMVDrGZ7WX2UPFxARByA4i7ELWqXHZxiXU2QqVFAQhhAi7kHVMJsWUMltK25Y7LFLJURCyjAi7kBMqHFbsYwi21Wyitsw+QSMShMmDCLuQM2rLk7tklFJMr3RgEt+6IGQdEXYhZ9gtZk6rdibMJq0tt0mWqSDkCKm4JOQUq9lEQ7WTzgEvQUPjtJpxhH8EQcgNIuxCzlFKMa3Cke9hCMKkQd6FBUEQSgwRdkEQhBJDhF0QBKHEEGEXBEEoMUTYBUEQSgwRdkEQhBJDhF0QBKHEEGEXBEEoMUTYBUEQSgyltZ74kyrVCRyc8BOPj1qgK9+DyAKlch0g11KIlMp1QHFcy0ytdd1YG+VF2IsBpdQmrfXSfI8jU0rlOkCupRApleuA0roWccUIgiCUGCLsgiAIJYYIe3J+nu8BZIlSuQ6QaylESuU6oISuRXzsgiAIJYZY7IIgCCWGCLsgCEKJIcKeBKXUD5RSu5RSW5VSjyqlqvM9pvGilLpSKbVbKfWOUuqr+R5PuiilmpRSLyqldiiltiul/ibfY8oEpZRZKdWqlHoy32PJBKVUtVLqt+H7ZKdS6t35HlO6KKX+Nvzd2qaU+rVSqqhbfomwJ2cdsEhrvRjYA9yZ5/GMC6WUGfgJ8AFgIfAxpdTC/I4qbQLAF7XWC4HlwF8V8bUA/A2wM9+DyAJ3A89orc8EllCk16SUagA+DyzVWi8CzMBN+R1VZoiwJ0Fr/ZzWOhD+73qgMZ/jSYNlwDta631aax/wAHBtnseUFlrrI1rrzeHf+wkJSEN+R5UeSqlG4CrgF/keSyYopaqAS4D/BtBa+7TWPfkdVUZYAKdSygK4gMN5Hk9GiLCnxu3A0/kexDhpAA5F/b+dIhXDaJRSs4AWYEN+R5I2PwK+DBj5HkiGzAY6gXvDbqVfKKXK8j2odNBadwD/ArQBR4BerfVz+R1VZkxqYVdKPR/2qY38uTZqm7sIuQLuz99IBQClVDnwMPAFrXVfvsczXpRSVwPHtdZv5HssWcACnAv8TGvdAgwCRTmPo5SqIfQ2Oxs4DShTSn0iv6PKDEu+B5BPtNaXj7ZeKXUrcDXwXl18Af8dQFPU/xvDy4oSpZSVkKjfr7V+JN/jSZMLgZVKqQ8CDqBSKXWf1roYRaQdaNdaR96cfkuRCjtwObBfa90JoJR6BLgAuC+vo8qASW2xj4ZS6kpCr8wrtdbufI8nDV4H5imlZiulbIQmgx7P85jSQimlCPlyd2qtf5jv8aSL1vpOrXWj1noWoc/jD0Uq6mitjwKHlFLzw4veC+zI45AyoQ1YrpRyhb9r76VIJ4IjTGqLfQx+DNiBdaHPmvVa69X5HVLqaK0DSqk7gGcJzfL/j9Z6e56HlS4XAjcDbyml3gwv+5rW+qk8jkmAvwbuDxsO+4Db8jyetNBab1BK/RbYTMjt2kqRlxeQkgKCIAglhrhiBEEQSgwRdkEQhBJDhF0QBKHEEGEXBEEoMUTYBUEQSgwRdkEQhBJDhF0QBKHE+P8w4nZ3jw2aowAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc3a5aa2eb8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f = gpflow.models.GPR(X, Y, gpflow.kernels.RBF(1))\n",
    "f.compile()\n",
    "gpflow.train.ScipyOptimizer().minimize(f, maxiter=notebook_niter(1000))\n",
    "full_lml = plot_model(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Upper bounds for sparse variational models\n",
    "As a first investigation, we compute the upper bound for models trained using the sparse variational GP approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T12:31:49.203722Z",
     "start_time": "2018-06-19T12:31:18.932370Z"
    },
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 "
     ]
    }
   ],
   "source": [
    "Ms = np.arange(4, notebook_niter(20, test_n=6), 1)\n",
    "vfe_lml = []\n",
    "vupper_lml = []\n",
    "vfe_hyps = []\n",
    "for M in Ms:\n",
    "    Zinit = X[:M, :].copy()\n",
    "    vfe = gpflow.models.SGPR(X, Y, gpflow.kernels.RBF(1), Zinit)\n",
    "    vfe.compile()\n",
    "    gpflow.train.ScipyOptimizer().minimize(vfe, disp=False, maxiter=notebook_niter(1000))\n",
    "    \n",
    "    vfe_lml.append(vfe.compute_log_likelihood())\n",
    "    vupper_lml.append(vfe.compute_upper_bound())\n",
    "    vfe_hyps.append({p.pathname:p.read_value() for p in vfe.trainable_parameters})\n",
    "    print(\"%i\" % M, end=\" \")\n",
    "vfe_hyps = pd.DataFrame(vfe_hyps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T12:31:49.334283Z",
     "start_time": "2018-06-19T12:31:49.205321Z"
    },
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'LML bounds for models trained with SGPR')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xd8VfX5wPHPkwGBhD1lhilLloCAgKDUUS24RxW3gFJr1dbWUaWtVltx1J8bRXFUodY9ELEyZIgQRVnKHgkz7ISEjOf3x/ckXEJyczNuzk3yvF+vS+4959xznntJznO+43y/oqoYY4wxRYnyOwBjjDGRzRKFMcaYoCxRGGOMCcoShTHGmKAsURhjjAnKEoUxxpigLFFUYSIyUUTeiIA4EkVERSSmhO8TEXlFRPaKyOJwxVdeRGSjiIwMYbtSfR/lQUSGishPYdr3qyLyYBn3cY+IvBRk/bUi8nVZjmFKzhJFBSvqZCIiw72Tx3sFlvfyls8OWKYi0rECwvXbEOAXQCtVHeB3MH4rj8SvqvNU9cTyiqm8qerfVfVGKJ+EKiJDRGSBiOwXkT0iMl9E+gesP0FEJotIiogcEpH1XsLrUiCGQ95jo4j8KeD9KiJp3rpkEXlcRKLL8h1EIksUkWUXMEhEGgUsuwb42ad4/NYW2KiqaSV9ox9X637zSmD2N+0RkbrAx8D/AQ2BlsBfgExvfSNgAVAbGArUAfoCc3AXKIHqq2oCcAVwv4icHbCul7fuNOAy4PpwfSa/2C9VZDkCvA9cDuBdmVwGvFmGfcaJyDQROSgiSSLSK2+FiHQVkdkisk9EVojIqIB1s0XkxoDXxxT5vSup8SKyxnv/MyIieXGLyCQR2S0i64FzAwPy9rXei2mDiFxZMGgRuQF4CZc4D4nIX7zlN4nIWu/q8EMRaVEgpgkisgZYU8g+864OrxORLV6V1ngR6S8iP3if4+mA7aNE5D4R2SQiO0XkNRGpF7B+jLcuVUTuLXCsKBH5k4is89ZPF5GGhf0Hhfh9nA3cA1zmfR/LAv6fHhKR+UA60N77fKu8/a0XkXEB+xkuIlsDXm8Ukd97n3+/97sSF7D+PBH53vtuFohIz4B1fbzfqYMiMg3If18h8W8SkZO951d6/w/dvdc3iMj73vPAUtNc7+c+7zMPCtjfJO//b4OInFPEYTsDqOpbqpqjqodVdaaq/uCtvx04AIxR1XXq7FPVV1T1/wrboaouBFYAPQpZtxaYD/Qu6nuotFTVHhX4ADYCIwtZPhzYCgwGvvGW/RL4HLgRmB2wrQIdQzjWRCALuBiIBX4PbPCexwJrcSefGsDpwEHgRO+9s4EbA/Z1LfB1gRg+BuoDbXClobO9deOB1UBr3JXcV972MUA87o8z7zgnAN2LiL/gMU8HduOu+mrirhTnFojpC++YtQrZX6K3zfO4k9qZQAYuOTfFXXHuBE7ztr/e+47aAwnAu8Dr3rpuwCFgmBfL40B23v8tcBuwCGjlrX8BeKtAHCX9PiYCbxRYNhvYDHT39heLS8wdAMFd5aYDfQN/zwr8Pi4GWnjf2ypgvLeuj/d9nAJE40q3G73PUwPYhDvZxuJ+x7KAB4uI/TXgTu/5i8A64OaAdbcX/IyB31OB34ks4CYvppuBFEAKOWZdIBWYCpwDNCiwfhEwsZi/ocD/KwFO9b7PMwr+LQJdgG15n6UqPaxEEWFUdQHQUEROBK7G/RGVxVJVfUdVs3AnszhgoPdIAB5R1SOq+j/cif+KEuz7EXVXYJtxySDvSupS4ElV3aKqe4CHC7wvF+ghIrVUdZuqrgjxeFcCU1Q1SVUzgbtxJY7EgG0eVtU9qno4yH7+pqoZqjoTSMOdwHeqajIwD3eCzDve46q6XlUPece7XFy11sXAx6o614vlz97nyjMeuFdVt3rrJwIXS+FVYqX9PvK8qqorVDVbVbNU9RM9eoU8B5iJq1opylOqmuL9X33E0f/HscALqvqNuivyqbhqm7zfn1jc/3OWqr4DfBvkGHNwSQsvlocDXp/mrQ/VJlWdrKo5uCRwAtCs4EaqegDXzqXAZGCXVwrN27YxsD1vexEZ5ZWcDorIzAK72w3swZVy/6SqXwasSxKRNFySnQ08W4LPUilYoohMrwO/AUYA7xWzbXG25D1R1VxcqaWF99jiLcuzCXdVHartAc/TcYmHvH0X2G9eDGm46rTxwDYR+US8hsMQtCiwr0O4K8bAmLcUfFMhdgQ8P1zI68DPsSlg3SbclWUzCnxG73OlBmzbFnjPO/Hsw51EcihwQivj95HnmM8sIueIyCKvem4frmTaOMj7i/p/bAvcmfcZvH215ujvT7J6l9KewO+qoDnAUBE5AVcSmA6c6iX5esD3xXzGQuNV1XTvaUJhG6rqKlW9VlVb4aqLWgBPeqtTcUkmb9sPVbU+rpRUo8CuGqtqA1XtqqpPFVjX1zv+ZbjSV3wJPkulYIkiMr0O3AJ8GvCHUFqt856Ia+hshSuqpwCt5djGzzZAsvc8DdfIl6d5CY65LfC43n7zqernqvoL3B/patzVXihScCcvAEQkHmgUEDO4q8fycszxcJ8jG5dYjvmMIlLbiyXPFuAcVa0f8IjzSi3HKMH3UdRny18uIjWB/wKTgGbeie9TXLVJSW0BHirwGWqr6lu4z99SRAL326bw3eTX36cDt+KqCw/gTvhjcdWLuYW9rRQxF0lVVwOvcrR94UvgfCmHDgBe6W06sBC4v6z7izSWKPwRKyJxAY9jqiNUdQOuOH5v4W8HoEaBfRTVJe9kEbnQO8bvcFUHi4BvcH+4d4lIrIgMB34FvO2973vgQhGpLa4r7g0l+HzTgd+KSCsRaQAEdidsJiKjvZN8Jq6ev7CTRGHeAq4Tkd7eCfHvuPacjSWIrSTeAm4XkXYikuAdb5qqZgPvAOeJ635ZA/grx/49PQ88JCJtAUSkiYiMLniAEn4fO4DEYk5sNXBtCLuAbK+h98wSfOZAk4HxInKKOPEicq6I1MGdELNx/8+xInIhUFwX5jm4knJeNdPsAq8L2oX7LtqXJngR6SIid4pIK+91a1zV6iJvk8eBBsDrItLB+4x1KFtj9CPATSJSkguriGeJwh+f4qo48h4TC26gql+rakqQfawosI/ritjuA1yReC8wBrjQq1M+gksM5+DqX58FrvauugCewPXC2oGrBy5Jz6vJuEb4ZUASrhE4TxRwB+5qfQ8uId4cyk5VdRauLeC/uCvaDng9xMJkCq50NxfXCSADd0WM144wAfi3F8teXLVenn8BHwIzReQg7uR0SiHHKMn38R/vZ6qIJBW2gaoeBH6LS9Z7gV97cZSYqi7BNRo/7e1rLa4xGe/350Lv9R7c79i7he0nwBxcF9S5RbwuePx04CFgvlf1NbCEH+Eg7jv/xmtDWAQsB+709r8b19aSAXztbf+9F1NIv5OFxPwj7vP8oTTvj1RybBWjMcYYcywrURhjjAnKEoUxxpigLFEYY4wJyhKFMcaYoKrEwGmNGzfWxMREv8MwxphKZenSpbtVtUlx21WJRJGYmMiSJUv8DsMYYyoVEQl2N30+q3oyxhgTlCUKY4wxQVmiMMYYE1SVaKMwxphgsrKy2Lp1KxkZGX6H4ou4uDhatWpFbGxsqd5vicIYU+Vt3bqVOnXqkJiYyLED3lZ9qkpqaipbt26lXbt2pdqHVT0ZY6q8jIwMGjVqVO2SBICI0KhRozKVpixRGGOqheqYJPKU9bNb1VN1s2UxrPsfSDREx0J0jYCfec+LWu49j4o5dnlsLfcwxlRJliiqg9xc+OlTWPAUbPkmPMdoPxz63wSdz4Zo+7UypqCEhAQOHTrkdxilYn/RVVnWYVj2Fix4Gvasg/pt4Jx/Qu8rIaYm5ByBnCzvcSTgtfc8N7vw5ce8JwvSdsKyt2HalVC3FfS7FvpeAwlN/f4GjKl2cnJyiI4uasLL0rFEURWlpcK3L8HiFyF9N7ToAxe/Al1HHXu1H126rnKFGn4P/PyZO+7/HoTZ/4Buo6D/jdBmEFTj+mFjAqkqd911F5999hkiwn333cdll13GhAkTOOussxg1ahQXXHABDRo0YMqUKUyZMoV169bx0EMP8cYbb/DUU09x5MgRTjnlFJ599lmio6NJSEhg3LhxzJo1i2eeeYYhQ4aUa8yWKKqSPeth4TPw3ZuQfRg6nQWn/hbanhr+E3V0DHT9lXvsXgNLprg4lv8XmnaH/jdAz8ugZkJ44zCmGH/5aAUrUw6U6z67tajLA7/qHtK27777Lt9//z3Lli1j9+7d9O/fn2HDhjF06FDmzZvHqFGjSE5OZtu2bQDMmzePyy+/nFWrVjFt2jTmz59PbGwst9xyC2+++SZXX301aWlpnHLKKTz22GPl+rnyWKKoCrYugfn/glUfuVJCz0th0K3QtIs/8TTuBGc/DKffBz++A99Ohk/ugC8egN5XQL8b/IvNGJ99/fXXXHHFFURHR9OsWTNOO+00vv32W4YOHcqTTz7JypUr6datG3v37mXbtm0sXLiQp556iqlTp7J06VL69+8PwOHDh2na1FXvRkdHc9FFF4UtZksUlVVuLvw8wzVQb14IcfVgyO1wyjio09zv6Jwa8XDyNdD3atj6rauWWvqqqxJLHOqqpbqcW75VYMYUI9Qr/4rWsmVL9u3bx4wZMxg2bBh79uxh+vTpJCQkUKdOHVSVa665hocffvi498bFxZV7u0Qgu4+issnKcCfbZwbA21fA/mQ4+xG4fSWMfCBykkQgEWg9AC580cV5xgOwdxP85xp48iSY/Qgc2OZ3lMZUiKFDhzJt2jRycnLYtWsXc+fOZcCAAQAMHDiQJ598Mr8qatKkSQwdOhSAM844g3feeYedO3cCsGfPHjZtCmmU8DKzEkVlkb4Hvn0ZFr8AabvghF5w0cvQ7fzK1R01oQkMvQNOvQ3WzHSljNkPw9xHoct5rpSROMQav02VdcEFF7Bw4UJ69eqFiPDPf/6T5s3dBd7QoUOZOXMmHTt2pG3btuzZsyc/UXTr1o0HH3yQM888k9zcXGJjY3nmmWdo27Zt2GMWVQ37QcKtX79+WiUnLsrNhZ0rIWkqfPcGZKVDx1+4BurEoVXnZJq6zmv8fgMy9kHjE+GMP7uGcWPKwapVq+jatavfYfiqsO9ARJaqar/i3uvLpaiIXAJMBLoCA1R1ScC6u4EbgBzgt6r6uR8x+iI3B7b/AJsWHH0c3gNReQ3Uv4Fm3fyOsvw16gBnPeQav5f/FxY+C9OucqWO0++vXCUmY6ogv/4ClwMXAi8ELhSRbsDlQHegBTBLRDqrak7Fh1gBso9Aynewab5LCpsXwZGDbl2DdnDiL6HtYOh4RmS2PZS32FrQ5yo46RKY8SfXkyvle3cPSHwjv6MzptryJVGo6ioodKCq0cDbqpoJbBCRtcAAYGHFRhgmWYddV9ZN891jy7fufgeAJl2g5yXunoe2g6FuC39j9VNMTTjvCWjRFz65E148DS573d04aIypcJFWpm8JLAp4vdVbdhwRGQuMBWjTpk34IyuNzINubKVNC2DjfEheCrlZgEDzk+Dka11SaDsY4hv7HW3k6TsGmnWH6VfDy2e55NHnSr+jMqbaCVuiEJFZQGH1Jfeq6gdl3b+qvgi8CK4xu6z7Kzcb5rn7GzYtgG3LQHPcSK0t+sCgW1yJofUpUKu+35FWDi37wtjZ8M718MEtLtme/QjE1PA7MmOqjbAlClUdWYq3JQOtA1638pZVDul74LVRrvG5VT8YeqcrLbTqb0NXlEV8Y7jqXfjfX127xfYf4NLXqnf1nDEVKNKqnj4E/i0ij+MaszsBi/0NqQSSk0Bz4ap3oN0wv6OpWqJj4Bd/dSWz9yfAC6fBpVNdIjbGhJUvd2aLyAUishUYBHwiIp8DqOoKYDqwEpgBTKhUPZ5SkgCxRtdw6n4B3PQl1KwDU38Fi56HKnAvkDHlSVXJzc0tt/35kihU9T1VbaWqNVW1maqeFbDuIVXtoKonqupnfsRXaslJ0LizO4mZ8GnaFcZ+BZ3OhBl/hHfHwpF0v6MyJqiNGzfSo0eP/NeTJk1i4sSJDB8+nNtuu43evXvTo0cPFi92lSgTJ05kzJgxDBo0iE6dOjF58uT89z766KP079+fnj178sADD+Tv/8QTT+Tqq6+mR48ebNmypdxij7Sqp8pL1ZUo2o/wO5LqIa4eXPYmzHsMvnoIdq5yXWgbtvM7MhPpPvsTbP+xfPfZ/CQ455FSvz09PZ3vv/+euXPncv3117N8+XIAfvjhBxYtWkRaWhp9+vTh3HPPZfny5axZs4bFixejqowaNYq5c+fSpk0b1qxZw9SpUxk4cGB5fTLABgUsPwdS4NAO10vHVIyoKDjtD3Dlf2D/ZnhxOKyZ5XdUxpTYFVdcAcCwYcM4cOAA+/btA2D06NHUqlWLxo0bM2LECBYvXszMmTOZOXMmffr0oW/fvqxevZo1a9YA0LZt23JPEmAlivKTkuR+trBEUeE6/cJ1oZ02Bt68GE6/F4bc6RKJMQWV4cq/LGJiYo5pN8jIyMh/XvDm47zXhS1XVe6++27GjRt3zLqNGzcSHx9f3mEDVqIoP8lJEBXjiqCm4jVsDzd8ASdd7KZinXYVZOz3Oypj8jVr1oydO3eSmppKZmYmH3/8cf66adOmAW5So3r16lGvXj0APvjgAzIyMkhNTWX27Nn079+fs846iylTpnDo0CEAkpOT84ceDxcrUZSXlCRo2g1i4/yOpPqqURsunAwt+8Hn98Dk0107hs2mZyJAbGws999/PwMGDKBly5Z06XL09zIuLo4+ffqQlZXFlClT8pf37NmTESNGsHv3bv785z/TokULWrRowapVqxg0aBAACQkJvPHGG2GduMiGGS8PqvCPtq7r5q/+5V8c5qiN893ESEfS4fxnofv5fkdkfBTJw4wPHz6cSZMm0a/fsaN9T5w4kYSEBH7/+9+Xy3HKMsy4VT2Vhz3rXTWHtU9EjsRTYdxcNyz7f66BmfdBdqbfURlTKVnVU3lI9hqyrcdTZKnbAq79BGbcDQv+D9Z+Cec/By16+x2ZMflmz55d6PKJEydWaBzBWImiPKQkQUwtaBKZRdtqLaYmnPc4/Hq6G4tr8unw1d/dXCDGmJBYoigPyUlwQk+biS2SdT4LJixykyLN+Qe8dDpsX+53VMZUCpYoyion2w0nbu0Tka9WA7jwBbj833Bwu7tBb+6j7v/QGFMkSxRltWu1m6XO2icqjy7nwi3fQLdR7p6Ll85wQ4AYYwpliaKsUr5zP61EUbnEN4KLp8AlU2H/FnhhGHz9BORWnsGKTeXy1FNP0bVrV668suhZGhMS3Lw1BQcQ9JtVqpdVShLUrOfuDDaVT/fz3ayDn9wOsybC6k9cz6jGnfyOzFQxzz77LLNmzaJVq1Z+h1JiVqIoq+Qk193SxhWqvBKawKWvw0UvQ+paeH4ILHzGShem3IwfP57169dzzjnnUK9ePSZNmpS/rkePHmzcuNG/4EJgJYqyyM6EHStg8G/8jsSUlYgbJypxKHz8OzcEyKqPYPQz0KiD39GZcvTzjoMczMgq133WiYulc7Oi56F5/vnnmTFjBl999RVPP/10uR67IthlcFlsXw65WdY+UZXUaeZ6RV3wAuxcCc+dCt+8AOU4W5gxlY2VKMoixe7IrpJEoNflbt7zD38Ln93llS6ehgaJfkdnyijYlX9FCDbceKSyEkVZJCdBfFOo29LvSEw41G3hJkUa9bS7V+bZwfDtyzZHtymTxMREkpLcRWZSUhIbNmzwOaLiWaIoi5QkV5ooMLmIqUJEoO8YuHkBtB4An9wBr58P+8pvPmJTvVx00UXs2bOH7t278/TTT9O5c2e/QyqWVT2VVuZB2PUTdL/Q70hMRajfGsa8B0tfdSPRPjcYznvCNYAbE4LAnk0zZ84sdJu8yYgSExPz582OBFaiKK1tywCFFn38jsRUFBHod50rXTTtBv+9AT681c15YUwVZomitGxo8eqrQVs3fPnQOyHpdTci7c7VfkdlTNhYoiitlCSo1wbiG/sdifFDdAyccT+MeRfSd7sBBpNet4buCFYVZvMsrbJ+dksUpZWcBC2t2qna63A6jJ/vGro//A28e5NrvzIRJS4ujtTU1GqZLFSV1NRU4uLiSr0Pa8wujbRU2LcJ+l3vdyQmEtRp5hq6v37cTYqUnASXvAIn9PI7MuNp1aoVW7duZdeuXX6H4ou4uLgyjTFliaI08kaMtfYJkycqGob9wQ0w+M4N8NJIOPMhGHCTdZ+OALGxsbRr187vMCotq3oqjZQkQOAEm3vZFNB2MIz/GtqPgM/+ANOugsN7/Y7KmDKxRFEayUluGOq4un5HYiJRfCP49TRXovj5c3h+GGxZ7HdUxpSaJYqSUnUlChsI0AQj4kYVvv5z93zK2fD1kza4oKmULFGU1IEUOLTD2idMaFqdDOPnQdfzYNYD8O9LIG2331EZUyKWKEoqb8RYK1GYUMXVc1Ounvs4bJjnhi7fMM/vqIwJmSWKkkpOgqgYaH6S35GYykQE+t8AN30JNevAa6Pgq4dtFj1TKViiKKmUJDfOT2zpb14x1Vjzk2DsbOh5Gcx5BKaOctWZxkQwSxQloeruobD2CVMWNRPggufh/OfchcfzQ2DNF35HZUyRLFGUxJ71kLHf2idM+ej9axg7B+qcAG9eDJ/eBUfS/I7KmONYoigJuyPblLcmneHGWTBgHCx+AZ4dBOtn+x2VMcewRFESyUkQUwuadPU7ElOVxNaCX/4TrvvMdZR4bbSbqztjv9+RGQP4lChE5FERWS0iP4jIeyJSP2Dd3SKyVkR+EpGz/IivSClJcEJPN8S0MeWt7WC4eT4M/i189zo8M9Dd2W2Mz/wqUXwB9FDVnsDPwN0AItINuBzoDpwNPCsi0T7FeKycbDernbVPmHCKrQVn/s1VR8XVg39fCu+Og/Q9fkdmqjFfEoWqzlTVbO/lIiBv/NvRwNuqmqmqG4C1wAA/YjzO7p8gK93aJ0zFaHkyjJsDp/0Rlr8DzwyAlR/4HZWppiKhjeJ64DPveUtgS8C6rd6y44jIWBFZIiJLKmSM+WS7I9tUsJiaMOIed99F3RYw/WqYNgYO7fQ7MlPNhC1RiMgsEVleyGN0wDb3AtnAmyXdv6q+qKr9VLVfkyZNyjP0wqUkQc160LB9+I9lTKDmJ8GN/4MzHnBtFs8MgGXTbNpVU2HC1iqrqiODrReRa4HzgDP06PyEyUDrgM1aecv8l5wELXpBVCQUwky1Ex0DQ++ALufBBxPgvbGw/L9w3hNQr9BCtzHlxq9eT2cDdwGjVDU9YNWHwOUiUlNE2gGdAP8H8s/OhB0rrNrJ+K9JZ7h+Bpz9CGyYC88OhKWvWunChJVfl8dPA3WAL0TkexF5HkBVVwDTgZXADGCCqvo/atr25ZCbZQ3ZJjJERcPAm+GWBW5e7o9uc/de7N3od2SmivLlhgBV7Rhk3UPAQxUYTvFsaHETiRq2h6s/hKSpMPPP7q7ukROh/01WRWrKlf02hSI5CeKbQL1WxW9rTEWKioJ+18GERdD2VPjsLnjlHNi9xu/ITBVitxiHIm/qUxG/IzGmcPVawZX/gR+mwWd/dJMj9bzUDThYu5H3aBjwvBHUqO131KaSCDlRiEjtAg3P1UPmQdj1E3S/wO9IjAlOBHpdDu1HwIw/wU+fend0F9HQHVOr8ARS1LJa9SEq1rWR2EVTtVJsohCRwcBLQALQRkR6AeNU9ZZwBxcRti0D1NonTOVRpxlc8op7npvjBhdMTy3isefo832b3M9QBiOMioXoWDeIYVRMEc9jXbfeQp/HuhsKO58N3c6HmBrh/U5MmYRSongCOAvXdRVVXSYiw8IaVSTJuyPbejyZyigq2isdNMT1Ng9BTtaxCSTvkbHPJZ6cLNcLMDfbjYFW6HPvZ+DznCw4ku5tk+P2ufy/8Pm90O9619ZSp3lYvw5TOiFVPanqFjm2qOl/l9WKkpIE9dpAfGO/IzGmYkTHulJJnWbhPU5uLqz7n5uHY84jMG+SK12cMg5a9bfqrQgSSqLY4lU/qYjEArcBq8IbVgRJToKWffyOwpiqJyoKOo10j9R1sHgyfP+mGwTxhN4uYXS/0OanjwChdI8dD0zADc6XDPQGqkf7RJpXb2vtE8aEV6MOcM4jcMcq+OUkyDoM798MT3SHL/8G+yNjJJ/qKpREcaKqXqmqzVS1qapeBVSPKd5s6lNjKlbNBBhwE0z4Bsa8D60HwLzH4MmTYPo1sGmBDVfig1Cqnv4PKHimLGxZ1ZOSBIgrBhtjKo4IdBjhHns3wrcvQdJrsPJ9N5rugLFw0iVuoicTdkUmChEZBAwGmojIHQGr6gKRMetcuCUnQeNOEFfX70iMqb4aJMKZD8Lwu+GH6bD4RfjwVvjifuh7DfS/Eeq3LnY3pvSCVT3VwN07EYMbwC/vcQC4OPyhRYCU76x9wphIUSPedaG9eQFc8zEkDoEFT8G/esLbV7rRdHNzrWoqDIosUajqHGCOiLyqqpsqMKYSSz+SzdJN5TyncNpuOFAH4gZCee/bVFo5uZCRlcPhrBxycpXcXCVX3SMnl/znublKjuLWo/nbKpCT45bl5kJO/vbuvYqiCrnqzneK9z516yHvGEfX56079nXo50st6s7tcNPAe8Y1P17l6GfLX8DR1259LeB3xLceQ7f9c+myej61VgZWfEAuUSgCCJr38LrcKlEgeMujjlmuSMA69/5jBe+2e/y3KQXWF3h9XDdgyf+Zv6/8+I7fx6E2Z5B42lVBYyqrUNoo0kXkUaA7kN9PTVVPD1tUkWDXz+5n4xP9jcOUC9WjJ/j0IzmkZ2aT7j0/fCSHtMxsDme55+n5j4BtMt17D2dF/i1EUeLOK+Kd8ELh1x0LEvhvQLgi4p7LsdvmnVPFe4MAsxhGDQbRX5bTmL2IgpAbcJrXow8NOP1r7vHrvddR5HrHOfa0f/wpvWBaCJ5087YvmAqO7vf4VHD8Nsdum5MZS2LQo5ZdKIniTWAabja68cA1QAVMUh262jViOLltw/Ld6dokiNlDt45lAAAa+UlEQVQAffpaP+5KIjdXSdl/mPW70li/6xAbdqexfnca63elsf1ABjm5xV85164RTZ24GBJqxlAnLpbmdeO85zEk1IylTpx7Hl8zhpgoIdp7RIkQEyVERQnRIkRHez8LbBMd5W1XYF20CFFRECXiPdzJMkqOLpP89e6nyPHbV2/n+R1AlRVKomikqi+LyG0B1VHfhjsw3yUnQdNuliQi0P7DWazfdYj1u9K8ZHD0eWZ2bv52CTVjaN8knn6JDWjVoBZ14mLzk0DduFgS4mKOJoWascTXjCYm2kbeN6agUBJFlvdzm4icC6QA5Xz5HmFUXUN29/P9jqTaOpKdy+Y96S4h7E5jw66jCSE17Uj+dtFRQpuGtWnfOJ6hnRrTrnEC7ZvE075JPE0SatpVtjHlIJRE8aCI1APuxN0/URe4PaxR+W3PejcAWgsbuqOiqCortx3g8xU7+GLlDn7ecfCYqqLGCTVp3zieX3RrRvsm8fkJoU3D2sRaKcCYsCo2Uajqx97T/cCI8IYTIfLuyLausWGVk6skbd7L58u38/nK7WzZc5gogX6JDblleIf8hNCucTz1asX6Ha4x1VYo81G0A24FEgO3V9VR4QvLZ8lJEBMHTavHSCUV6Uh2LgvW7c4vOew+lEmN6ChO7diICcM7MrJbMxon1PQ7TGNMgFCqnt4HXgY+AnKL2bZqSEmC5j3dcMumzNIys5nz8y4+X7Gd/63eycGMbOJrRDO8S1PO6t6cESc2oU6cfdfGRKpQEkWGqj4V9kgiRU62m9Wu79V+R1Kp7U07wqxVO/h8xQ7mrdlFZnYuDeNrcE6P5pzdozmDOzQmLrZ6jARjTGUXSqL4l4g8AMwEMvMWqmpS2KLy0+6fICvd2idKYdv+w8xcsYPPV2znmw17yMlVWtSL44oBbTi7R3P6tW1g3U+NqYRCSRQnAWOA0zla9aTe66rHpj4tkSPZuUxdsJGPf9zGsi37AOjYNIHxp7Xn7O4n0KNlXeuiakwlF0qiuARor6pHit2yKkhJgpp1oWEHvyOJeAcysrj5jaXMX5tKr1b1+MNZJ3JW9+Z0bJrgd2jGmHIUSqJYDtQHdoY5lsiQnAQtertpGk2Rtu/P4NpXFrN25yEmXdKLi09u5XdIxpgwCSVR1AdWe8N2BLZRVL3usdmZsGMFDJrgdyQR7aftB7n2lcUczMjmlev6M7RTE79DMsaEUSiJ4oGwRxEpti+H3Cxrnwhi4bpUxr6+hFqx0UwbN5DuLer5HZIxJsxCuTN7TkUEEhFSvIZs6/FUqA++T+YP//mBto1q8+r1A2hZ36ahNKY6CDYV6teqOkREDnLsIOsCqKpWvflBk5MgvgnUs/r2QKrKi3PX8/BnqxnQriGTx/SjXm27Qc6Y6iLYDHdDvJ91Ki4cn6UkudKEdefMl5Or/PWjFUxduInzep7AY5f2omaM3ShnTHVSbNceEXk9lGWVXuZB2PWTtU8EyMjK4eY3ljJ14SbGDmvPU5f3sSRhTDUUSmN298AXIhIDnByecHy0bRmg1j7h2ZN2hBunfst3W/bxwK+6cd2p7fwOyRjjk2BtFHcD9wC1RORA3mLgCPBiBcRWseyO7HybU9O55pXFpOw7zHNX9uXsHif4HZIxxkfB2igeBh4WkYdV9e4KjMkfKd9BvTYQ39jvSHy1bMs+bpj6Ldm5yps3nkK/xKo9maExpnih3H78sYjEA4jIVSLyuIi0DXNcFS8lCVpW7xnt/rd6B5e/uIi42Gj+e/NgSxLGGCC0RPEckC4ivXDToa4DXgtrVBUtfQ/s3Vit2yf+/c1mbpy6hI5NE3j3lsF0aGLjNRljnFASRbaqKjAaeFpVnwGqVpfZlOrbPqGqPDbzJ+5570eGdW7C22MH0rROnN9hGWMiSCi9ng56DdtjgKEiEgVUrbutkr05sk/o5W8cFexIdi5/evcH3k1K5vL+rXnw/B42X4Qx5jihnBUuww0GeL2qbgdaAY+W5aAi8jcR+UFEvheRmSLSwlsuIvKUiKz11lfMJX5KEjTqBHHVZ9yigxlZ3DD1W95NSuaOX3Tm4QtPsiRhjClUsWcGLzn8F8ib8X438F4Zj/uoqvZU1d7Ax8D93vJzgE7eYyyufST8kpOqVbXT9v0ZXPrCIhauS+XRi3vy2zM62eRCxpgihXJn9k3AO8AL3qKWwPtlOaiqHgh4Gc/RsaRGA6+pswioLyLh7cR/IAUOba82Ddkbd6dx4bPz2ZyaxsvX9ueSfq39DskYE+FCaaOYAAwAvgFQ1TUi0rSsBxaRh4Crgf3ACG9xS2BLwGZbvWXbCnn/WFypgzZt2pQ+kGp0o52qcve7P5J2JIdp4wbRo2X1qWozxpReKJXSmYHToHpDeGiQ7fO2myUiywt5jAZQ1XtVtTXwJvCbkgauqi+qaj9V7dekSRkmzklJgqgYaH5S6fdRScxcuYOF61P5/ZmdLUkYY0IWSolijojkDeXxC+AW4KPi3qSqI0OM4U3gU9wESclAYF1IK29Z+CQnQdOuEFu151bIzM7h75+uonOzBK4YUIYSmDGm2gmlRPEnYBfwIzAOd1K/rywHFZFOAS9HA6u95x8CV3u9nwYC+1X1uGqncqPqhu6oBu0TUxdsZFNqOn8+r5v1bjLGlEgoM9zlApO9R3l5REROBHKBTcB4b/mnwC+BtUA6cF05HvN4e9ZDxr4q3z6x+1Am//flWs7o0tTmtzbGlFgoVU/lTlUvKmK54hrPK0aKd6NdFS9RPDbzZw5n5XDPuV39DsUYUwlV7zqIDqfD5W+5NooqamXKAaZ9u5mrByXa+E3GmFIpVaIQkUnlHYgvajeELr+E6Ko1IkkeVeVvH6+kXq1YbjujU/FvMMaYQpS2RHFpuUZhwiKvO+wdv+hMvdpVMxkaY8KvtInCxnuIcHndYTs1te6wxpiyCTYValGz1giWKCJeXnfY164fYN1hjTFlEqzX01LcHdiFJYWs8IRjykNed9jTuzRlWGfrDmuMKZtgc2a3q8hATPl5/AvXHfZe6w5rjCkHwaqegt5coKpJ5R+OKatV2w7w9uLNXDu4nXWHNcaUi2BVT0uA5bj5J+DYKigFTg9XUKZ0VJW/frSSutYd1hhTjoIlijuAi4HDwNvAe6p6qEKiMqXyhdcd9q+ju1t3WGNMuSmyO4yqPqmqQ4BbcSO6fiki00Wkd4VFZ0KWmZ3DQ1532F9bd1hjTDkKZSrU9cAHwEzcBEadwx2UKTkbHdYYEy7BGrPbA5fjhgHfgqt++ruqHq6g2EyIrDusMSacgrVRrAV+wJUmDgBtgJtFXJu2qj4e9uhMSPK6w97zS+sOa4wpf8ESxV85OuVpwX6WxU6FaipGXnfYawYn0rGpdYc1xpS/YDfcTSxqnYj8LizRmBLJGx3WusMaY8KptK2ed5RrFKZUvli5gwXr3Oiw9WvX8DscY0wVZaPHVlLWHdYYU1FKmyisjcJnry3YxKbUdO6z7rDGmDAL1j32IIUnBAFqhS0iU6zUQ5k89eUaRpzYhNOsO6wxJsyCNWbXqchATOgeyx8dtpvfoRhjqgGrs6hk8rrDjhnU1rrDGmMqhCWKSkRVefAT6w5rjKlYligqkVmrdjJ/bSq3j7TusMaYimOJopLIzM7hoU9W0rFpAr8+xbrDGmMqjiWKSuK1BZvY6I0OG2vdYY0xFcjOOJWAdYc1xvjJEkUl8PgXP5Nu3WGNMT6xRBHhVm07wFuLNzNmoHWHNcb4wxJFBFu8YQ8T3kyibq1YfjfSusMaY/wRbD4K45PUQ5k8/Nlq3lm6lZb1a/Hsr/tad1hjjG8sUUSQ3Fxl+pItPDJjNYcysrl5eAduPb0jtWvYf5Mxxj92BooQq7Yd4L73l7N0014GtGvIQ+f3oFMzG27LGOM/SxQ+S8vM5slZPzNl/kbq1Ypl0iW9uKhvS/LmJjfGGL9ZovCJqvL5iu385aOVbNufwRUD2vDHs0+0tghjTMSxROGDLXvSuf+D5Xz10y66nlCXp3/dl5PbNvA7LGOMKZQligp0JDuXyfPW89SXa4iJEu47tyvXDk60GeqMMRHNEkUFWbBuN39+fznrdqVxTo/m3P+rbpxQzyYKNMZEPl8ThYjcCUwCmqjqbnEtuP8CfgmkA9eqapKfMZbVroOZ/P3TVbz3XTKtG9bilev6M+LEpn6HZYwxIfMtUYhIa+BMYHPA4nOATt7jFOA572elk5Or/HvxZh6dsZrDWTncenpHJozoSFxstN+hGWNMifhZongCuAv4IGDZaOA1VVVgkYjUF5ETVHWbLxGW0vLk/dz7/nKWbdnHoPaN+Nv5PWycJmNMpeVLohCR0UCyqi4rcL9AS2BLwOut3rLjEoWIjAXGArRpExkT+agqj3y2msnz1tMwvgZPXtab0b1b2D0RxphKLWyJQkRmAc0LWXUvcA+u2qnUVPVF4EWAfv36aVn2VV6+XrubF+au56K+rbj/vG7Uqx3rd0jGGFNmYUsUqjqysOUichLQDsgrTbQCkkRkAJAMtA7YvJW3rFJ45qu1NKtbk79f2IOaMdYWYYypGiq8A7+q/qiqTVU1UVUTcdVLfVV1O/AhcLU4A4H9laV9YummvSxav4ebhra3JGGMqVIi7T6KT3FdY9fiusde5284oXtu9lrq147ligGR0V5ijDHlxfdE4ZUq8p4rMMG/aEpn9fYDzFq1k9tHdia+pu9fqTHGlCsbO6IcPDd7HfE1orlmcFu/QzHGmHJniaKMNqWm8dGyFK4c2NZGfjXGVEmWKMrohbnriYmK4oYh7fwOxRhjwsISRRnsOJDBO0u2cnG/VjSrG+d3OMYYExaWKMrg5a83kJ2by/hhHfwOxRhjwsYSRSntSz/CG4s28ateLWjTqLbf4RhjTNhYoiilVxdsJP1IDjcPt9KEMaZqs0RRCmmZ2by6YCMjuzalS/O6fodjjDFhZYmiFN5avJl96VncMqKj36EYY0zYWaIooczsHCbPW8+g9o3o26aB3+EYY0zYWaIooXeTktlxIJNbRljbhDGmerBEUQLZObk8P2cdPVvVY0jHxn6HY4wxFcISRQl8unw7m1LTuWV4B5u1zhhTbViiCJGq8uxXa+nQJJ4zuxU2cZ8xxlRNlihC9NVPO1m9/SC3DO9IVJSVJowx1YclihCoKs98tY6W9WsxqncLv8MxxpgKZYkiBIs37GHppr2MO609sdH2lRljqhc764XgmdnraJxQg0v7tfY7FGOMqXCWKIrx49b9zP15F9cPaUdcbLTf4RhjTIWzRFGM5+aspU5cDFcNtGlOjTHVkyWKINbuPMRny7dz9aC21I2L9TscY4zxhSWKIF6Ys46aMVFcd6pNc2qMqb4sURQhed9h3vsumcv7t6FxQk2/wzHGGN9YoijC5LnrAbhpWHufIzHGGH9ZoijE7kOZvP3tZi7o05KW9Wv5HY4xxvjKEkUhXpm/gczsXMbbNKfGGGOJoqADGVm8tnAT5/RoTocmCX6HY4wxvrNEUcAbizZxMCObW4bbNKfGGAOWKI6RkZXDlK83MKxzE3q0rOd3OMYYExEsUQSYvmQLuw8dYYK1TRhjTD5LFJ6snFxemLOek9s2YEC7hn6HY4wxEcMShefD71NI3neYCSNsmlNjjAlkiQLIzVWem7OOLs3rMOLEpn6HY4wxEcUSBTBz5Q7W7jzELSM6WmnCGGMKqPaJQlV5dvZaEhvV5tyTTvA7HGOMiTjVPlF8vXY3P2zdz7jTOhAdZaUJY4wpqNonime/WkezujW5sG9Lv0MxxpiIVK0TRdLmvSxcn8pNQ9tTM8amOTXGmML4kihEZKKIJIvI997jlwHr7haRtSLyk4icFc44VGFY5yZcMaBNOA9jjDGVWoyPx35CVScFLhCRbsDlQHegBTBLRDqrak44Aji5bQNeu35AOHZtjDFVRqRVPY0G3lbVTFXdAKwF7ExujDE+8jNR/EZEfhCRKSLSwFvWEtgSsM1Wb5kxxhifhC1RiMgsEVleyGM08BzQAegNbAMeK8X+x4rIEhFZsmvXrnKO3hhjTJ6wtVGo6shQthORycDH3stkoHXA6lbessL2/yLwIkC/fv209JEaY4wJxq9eT4G3QF8ALPeefwhcLiI1RaQd0AlYXNHxGWOMOcqvXk//FJHegAIbgXEAqrpCRKYDK4FsYEK4ejwZY4wJjS+JQlXHBFn3EPBQBYZjjDEmiEjrHmuMMSbCiGrlbwcWkV3AplK+vTGwuxzDCQeLsewiPT6I/BgjPT6I/BgjLb62qtqkuI2qRKIoCxFZoqr9/I4jGIux7CI9Poj8GCM9Poj8GCM9vqJY1ZMxxpigLFEYY4wJyhKFd9NehLMYyy7S44PIjzHS44PIjzHS4ytUtW+jMMYYE5yVKIwxxgRlicIYY0xQ1T5RiEi0iHwnIh8Xv3XFE5H6IvKOiKwWkVUiMsjvmAKJyO0issIbGfgtEYmLgJimiMhOEVkesKyhiHwhImu8nw2C7cOH+B71/o9/EJH3RKS+X/EVFWPAujtFREWksR+xBcRRaIwicqv3Xa4QkX9GUnwi0ltEFnkzey4RkUox3061TxTAbcAqv4MI4l/ADFXtAvQigmIVkZbAb4F+qtoDiMbNUOi3V4GzCyz7E/ClqnYCvvRe++VVjo/vC6CHqvYEfgburuigCniV42NERFoDZwKbKzqgQrxKgRhFZARuArReqtodmFTI+yrKqxz/Hf4T+Iuq9gbu915HvGqdKESkFXAu8JLfsRRGROoBw4CXAVT1iKru8zeq48QAtUQkBqgNpPgcD6o6F9hTYPFoYKr3fCpwfoUGFaCw+FR1pqpmey8X4YbY900R3yHAE8BduAE9fVVEjDcDj6hqprfNzgoPzFNEfArU9Z7XIwL+XkJRrRMF8CTulz7X70CK0A7YBbziVY+9JCLxfgeVR1WTcVdsm3ETUO1X1Zn+RlWkZqq6zXu+HWjmZzDFuB74zO8gCvImHUtW1WV+xxJEZ2CoiHwjInNEpL/fARXwO+BREdmC+9vxu+QYkmqbKETkPGCnqi71O5YgYoC+wHOq2gdIw98qk2N49fyjcQmtBRAvIlf5G1Xx1PUJ9/2KuDAici9uiP03/Y4lkIjUBu7BVZdEshigITAQ+AMwXUTE35COcTNwu6q2Bm7Hqy2IdNU2UQCnAqNEZCPwNnC6iLzhb0jH2QpsVdVvvNfv4BJHpBgJbFDVXaqaBbwLDPY5pqLsyJswy/vpW5VEUUTkWuA84EqNvBucOuAuCJZ5fzOtgCQRae5rVMfbCryrzmJcbYGvje4FXIP7OwH4D2CN2ZFMVe9W1VaqmohrgP2fqkbU1bCqbge2iMiJ3qIzcJM6RYrNwEARqe1dtZ1BBDW2F/Ah7o8U7+cHPsZyHBE5G1cNOkpV0/2OpyBV/VFVm6pqovc3sxXo6/2ORpL3gREAItIZqEFkjdaaApzmPT8dWONjLCHza4Y7E7pbgTdFpAawHrjO53jyqeo3IvIOkISrLvmOCBiiQETeAoYDjUVkK/AA8AiuGuIG3JD0l0ZYfHcDNYEvvJqSRao6PpJiVNWIqiYp4nucAkzxuqQeAa7xq3RWRHw3Af/yOn9kAGP9iK2kbAgPY4wxQVXbqidjjDGhsURhjDEmKEsUxhhjgrJEYYwxJihLFMYYY4KyRGHCxhth9LGA178XkYnltO9XReTi8thXMce5xBu196sCy1t4XYNLsq9rReTpUsbxVxEZWZr3llUoxxaR4SISqTdbmjKy+yhMOGUCF4rIw6oaMTc9iUhMwAB8xbkBuElVvw5cqKopQNgTVcDxfBs6I8RjDwcOAQvCG43xg5UoTDhl427Au73gioIlAhE55P0c7g3m9oGIrBeRR0TkShFZLCI/ikiHgN2M9Mb0/9kbuytvfpFHReRbb26HcQH7nSciH1LI3e0icoW3/+Ui8g9v2f3AEOBlEXm0wPaJefMMeCWFd0Vkhrj5Lv4ZsN11XnyLccPGBP383vM/erEsE5FHCm4vIhtF5C8ikuRt18Vb3kTcXBsrvAEkN0khc0aIyCERecLb7ksRaeItz5srIW9OjAahHFtEEoHxwO3i5lkY6pXElnufYW7BGEzlYonChNszwJXihkwPVS/ciacrMAborKoDcMPB3xqwXSJurJxzgefFTZp0A24U2/5Af+AmEWnnbd8XuE1VOwceTERaAP/ADanQG+gvIuer6l+BJbixl/5QTMy9gcuAk4DLRKS1uDGl/oJLEEOAbsV9cBE5BzfQ4imq2oui5yvYrap9geeA33vLHsANRdMdNy5YmyLeGw8s8bab470P4DXgj96cGD8GLA96bFXdCDwPPKGqvVV1Hm7wwLO8zzCquM9tIpslChNWqnoAdwL6bQne9q2qbvPmFFgH5A1d/iMuOeSZrqq5qroGN7xJF9ykOleLyPfAN0AjoJO3/WJV3VDI8foDs73BDfNGbh1WgnjBTYq0X1UzcCWWtsApAfs9AkwLYT8jgVfyxntS1cLmhICjA8st5eh3MgQ3wCWqOgPYW8R7cwNieQMY4iXy+qo6x1s+laK/g8KOXdB84FURuQk3oZWpxCxRmIrwJO5KP3AujWy83z8RicIN3pYnM+B5bsDrXI5tVys4/owCAtzqXdn2VtV2AXNkpJXpUwQXGHMOxbf/Bfv8JTleKMcqTknH8Sn22N44VfcBrYGlItKo9OEZv1miMGHnXRVPxyWLPBuBk73no4DYUuz6EhGJ8tot2gM/AZ8DN4tILLgRRKX4yZ4WA6eJSGMRiQauwFXJlNU33n4befFcErBuI4V//i+A68TN/4CINCzB8ebjDXYoImcCRc0LHsXRhvhfA1+r6n5gr4gM9ZaPoWTfwUGgTt4LEemgqt94DeG7cAnDVFLW68lUlMeA3wS8ngx8ICLLgBmU7mp/M+4kXxcYr6oZIvISrjokSUQEd5IKOu2pqm4TkT8BX+FKJJ+oapmHIff2OxFYCOwDvg9YXejnV9UZItIbWCIiR4BPcRMGheIvwFsiMsY75nbcCbygNGCAiNyHm5fjMm/5Nbi2ntqUfKTij4B3xM2CdyuuYbsT7vv8EojkWfFMMWz0WGOqCBGpCeSoaraIDMLNjNi7kO0OqWpCxUdoKisrURhTdbTBzbkRhZuL4Saf4zFVhJUojDHGBGWN2cYYY4KyRGGMMSYoSxTGGGOCskRhjDEmKEsUxhhjgvp/3JWPosZI7uUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc3a92b81d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(Ms, vfe_lml, label=\"lower\")\n",
    "plt.plot(Ms, vupper_lml, label=\"upper\")\n",
    "plt.axhline(full_lml, label=\"full\", alpha=0.3)\n",
    "plt.xlabel(\"Number of inducing points\")\n",
    "plt.ylabel(\"LML estimate\")\n",
    "plt.legend()\n",
    "plt.title(\"LML bounds for models trained with SGPR\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We see that the lower bound increases as more inducing points are added. Note that the upper bound does _not_ monotonically decrease! This is because as we train the sparse model, we also get better estimates of the hyperparameters. The upper bound will be different for this different setting of the hyperparameters, and is sometimes looser. The upper bound also converges to the true lml slower than the lower bound."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Upper bounds for fixed hyperparameters\n",
    "Here, we train sparse models with the hyperparameters fixed to the optimal value found previously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T12:32:58.327730Z",
     "start_time": "2018-06-19T12:31:49.335336Z"
    },
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 "
     ]
    }
   ],
   "source": [
    "fMs = np.arange(3, notebook_niter(20, test_n=5), 1)\n",
    "fvfe_lml = []  # Fixed vfe lml\n",
    "fvupper_lml = []  # Fixed upper lml\n",
    "for M in fMs:\n",
    "    Zinit = vfe.feature.Z.read_value()[:M, :].copy()\n",
    "    Zinit = np.vstack((Zinit, X[np.random.permutation(len(X))[:(M - len(Zinit))], :].copy()))\n",
    "    init_params = {p.pathname:p.read_value() for p in vfe.parameters} # TODO: provide convenience function\n",
    "    init_params['SGPR/feature/Z'] = Zinit\n",
    "    vfe = gpflow.models.SGPR(X, Y, gpflow.kernels.RBF(1), Zinit)\n",
    "    vfe.assign(init_params)\n",
    "    vfe.kern.variance.set_trainable(False)\n",
    "    vfe.kern.lengthscales.set_trainable(False)\n",
    "    vfe.likelihood.set_trainable(False)\n",
    "    vfe.compile()\n",
    "    gpflow.train.ScipyOptimizer().minimize(vfe, disp=False)\n",
    "    \n",
    "    fvfe_lml.append(vfe.compute_log_likelihood())\n",
    "    fvupper_lml.append(vfe.compute_upper_bound())\n",
    "    print(\"%i\" % M, end=\" \")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T12:32:58.517509Z",
     "start_time": "2018-06-19T12:32:58.329700Z"
    },
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc3517f1e10>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEKCAYAAADTgGjXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xl8lPW5///XlQXCkrDvEILIIiAEBVwqiNa6HD3iUlutW2t/LrWtfnvantZjW2lP7eJyaltte7S11qOn6nHD2ta1BtSKCIjIYtgECQIJIGtIyHL9/rjvJENIYDLJ5J4k7+fjcT9m7s993zPXsMw1n+X+fMzdERERSURa1AGIiEjbpSQiIiIJUxIREZGEKYmIiEjClERERCRhSiIiIpIwJREREUmYkoiIiCRMSURERBKWEXUAyda3b1/Py8uLOgwRkTZl0aJF29y935HOa/dJJC8vj4ULF0YdhohIm2JmG+I5T81ZIiKSMCURERFJmJKIiIgkTElEREQSpiQiIiIJUxIREZGEKYmIiEjC2v19Igl7+7/Bq6H3UcHWczhkdIo6KhGRlKIk0piFf4SSlXX7lgY9hkLvkXWJpWbrlQeZWZGFKiISFSWRxtz4FpRuhx3rgm372rrny56Esl0xJ1uYYEbUSzAjgwTTqWtUn0JEJKmURBpjBt36BtuwaYceL90BOz6sSyw7wiSz8i9B8omVPThIJj1zodfw4LFnbtBEljME0vXXICJtk769EtW1d7ANPf7QY/t3wicxCWb7Oti5ATa8Ce8/EfS11LD0IJEckmBqksxgSEtvvc8lItIESiLJ0KUndJkMgycfeqzyAOzeBDs/ChLLzo/qtrWvwZ7NgNedn5ZRL8mEiaX7AOjWr+5RtRkRiYC+eVpbRqew72REw8cry2FX0aEJ5pMNsPoV2LulgYssqBV1HwDd+0O3/sFj9wHh1q/ueZfekKaR3SLSMpREUk1GZ+gzMtgaUlEWJJK9xbB3a/hYDPuK68p2vB08VpYder2lhzWYMNH0yoNB+TBoEvQ/BtIzk/rxRKR9URJpazKzgi/+XnmHP88dyvfAvpIw2WyFveHzmoSzZwt89Da88/vgmvTOMGA8DA6TyqB86D9O98eISKPaXBIxs7OBXwLpwO/d/WcRh5SazCArJ9gaq9UAVFcHnf+bl8DH78Lm9+D9p2Dhg8Hx9E5BIhk0KUwu+UGiyejcOp9DRFKaufuRz0oRZpYOrAI+AxQB7wCXufuKxq6ZMmWKa2XDJqquDkaXbV4SJJWPlwTPa+6NScsImr4G5YeJZXKQWHTDpUi7YWaL3H3Kkc5razWRacAad18HYGaPAbOARpOIJCAtra5fZsLFQZk7fLL+4MTywfPw7v8Exy0d+o0Jkkv/ceF2TDCaTB35Iu1WW0siQ4CNMftFwAkRxdKxmNWNKht/YVDmHowc2/xekFy2LoeN78Cyp+quy+wG/cfGJJdjoP/4oFPfLJrPIiItpq0lkbiY2XXAdQC5ubkRR9OOmQX3rvQaDuPOrysv2w0lhVC8AopXQvFyWPUivPtI3TldetcllQFhzaXf2OAeGxFpM9paEtkEDIvZHxqWHcTd7wfuh6BPpHVCk1pZOTBsarDF2lsSTGpZvDJIMFtXwHuPwYE9defkDAlrK8cESaXfWOg7OnhNEUk5bS2JvAOMMrMRBMnjUuAL0YYkceveL9hGzKgrcw9urqxJLDU1l/VvHHyfS86QIJn0Gxv0vdQ8du3d+p9DRGq1qSTi7pVm9jXgRYIhvg+6+/KIw5LmMIOew4Jt9Jl15dVVwV37JYVQ8kHd4+KHoWJf3Xnd+h+cVGpqL936qs9FpBW0qSG+idAQ33amuhp2F8Ukl5oEUwjlu+vO69KrLrEMmgRHzQym5xeRuLTXIb7S0aWl1c1yPOozdeXuwR34sbWWkkJYMQcWPRSc03N4kExqNjWFiTSbkoi0D2aQMyjYRp5WV+4eLCi27jVYVwDLn4HFfwIMBk2Eo04LEkruSbpZUiQBas6SjqWqMpjeZV1BkFg2LoDqCsjIgtwTw1rKaTBwom6SlA4t3uYsJRHp2Mr3woZ/1iWV4nDygy694ahT65JKr+ERBinS+tQnIhKPzt2DUWE1I8P2bIF1c+uSyvJngvJeI4JmsrzpQY0lZ3BkIYukEtVERBrjDttWBStOriuA9a/Dgb3BsR65MGxakFCGnRBMQKlljKUdUU1EpLnMwntPxsCJN0BVBWxZGqzBsvFt2PAmLHsyOLdTdxg6BYadGCSXoVN1l710CKqJiCSqZgLKjQtg4/wgsWxdDl4NlhZMNBlbW+mZqxsgpc1QTUQk2WInoJx4SVBWths2LQwSy0fzYekTsPAPwbHsQUFSGRYmlUETtRyxtHlKIiItKSsHRp4ebBBM31K8IkgoG8NmsBVzgmMZWTB4ctD0VdMElj0wuthFEqDmLJHWtntz2Pz1DhQtCNZjqToQHOuRG8x+PDRMKgOP1Rr3Egk1Z4mkqpxBwcJeNYt7VZbD5qVBQqlpBqtZ2CsjK1iGeNjUIKkMnRZcL5IiVBMRSUW7NkHRO8G2cUGwcmRtbWVYmFDCZrCBE1VbkRanmohIW9ZjSLCNvyDYryyHLe8HCaWmxrL86eBYemcYOCFcyOuYcDnicUFHvkaDSZKpJiLSVu3+uK6msmVpsKDXvpK645171K1vH5tcuvVTcpEjUk1EpL3LGQzjZgVbjX3bY5YgXhlMib9iDux/qO6cLr0PXoK4Jsl069PqH0HaPiURkfakWx/odgrknVJX5g57iw9NLkufOHghr279g9pKz1xIywzuYUnLhPSMmP2MmPL6+xmHXpfZLVi3pWsfyOqhGlA7lHJJxMzuBP4VOACsBb7k7jvNLA9YCRSGp8539xsiCVKkLTGD7AHBdtTMunL3oEmsNrl8EDxf82owxUt1RTB1fnVFsE8zm77TMoJaULe+QVLp2hu61jzvE5b3rtvv2ldrvLQBKZdEgJeBW8L11H8O3AJ8Jzy21t3zowtNpPncnapqp8od9+C73Amfh8eDR+Awxzw4WHcuYdlB+3XvGVteF0sv6HUy9DoZH9PwOXUnV0N1BVZViVVXQHXwaNVBoql9jD1eWUpa2Sek799O+v4dpO3fTnrZDtL27yB91/KgrOwTrJEEVZ3Rlaouvanq0ofqTjl4eqeDNtJi9tMy8fTOh5wTe4za55lA/VpRA7WkBmpODUYaVw0r3lqYxbxm7HPwevsHvWb9MjOGjDqOtPTkTgyacknE3V+K2Z0PfLY5r1d6oJJFG3Y0Lyhp1yqqnH3llewrr6T0QBX7DlRSWh487iuvZN+BqtpjZRVVQQIIt8owGVRVBc+rvaa8mupqguO1z6vD5BH1J25NaUCfcGuYUU139pNDKdmUkmOl5Ni+4HlFKdn7S8mhlG5WRiYVZFBOBpVkWhWZVAXPqSKDKjKsqtU+WVuw8Yrn6Nw5ubW5lEsi9VwDPB6zP8LM3gV2A99z99ejCUtS3e79lWzcUcrGT0rZtvcApQcq2VcvMZSGiePAEb7V0wy6dUqna+cMOmekk5FmpIdbRprROS2N9M51ZelpaWQYpKcbaZZGRjqkW9pB16SnGWnhL0ezcKt5HvPr0sxqy2v3rea3ZnhmzLXh7kFPrN4v4Nrfqw38erb6P3RbQP33r/9e8SgNtyPyKqy6kjSvwqorSKuuwLyS9LDWlOaVWHUl5gcnm8ZqQg28QQNFR742/tdv6DX9oDKLLW+oLOZxTEby52aLJImY2StAQ5ME3eruc8JzbgUqgUfDY5uBXHffbmbHA8+a2Xh3313/RczsOuA6gNzcXI4f3jsZH0NSwN7ySlZv3cOqrXso3LKXVeHz4j3ltedkpBk5XTLJycogp0sm/bOzyOmSQU5WZm15dlZmvbK6/a6d0hv8whWRiJKIu59xuONm9kXgPODTHjbmuns5UB4+X2Rma4HRwCE3gbj7/cD9ENwn0qLBSyTKKqpYU7yX1cV1yaJwyx427dxfe05WZhqjB2QzfVQ/xgzszugB2YwekM2gHllKAiJJknLNWWZ2NvDvwKnuXhpT3g/Y4e5VZnYUMApYF1GYkkTb9pYzf912Vm3dy6otQc1i/fZ9VIc/BzLTjZH9unP88F584YRcRvXvzpiB2Qzr1ZW0NCULkdaUckkEuBfoDLwc/nqsGco7A/iRmVUA1cAN7q4e83akvLKKB99Yz6//sZrSA1WkGeT17caYgdmcN2kwYwZkM2Zgd4b36UZmelrU4YoIKZhE3P3oRsqfAp5q5XCklbxWWMyP/rKCD7ft4zPjBvC1045mzMBssjK1brlIKku5JCIdy4bt+/jP51fwyspijurbjYe+NJWZY/pHHZaIxElJRCJReqCS+15bwwPzPiQz3bjlnLF86VMj6JShZiqRtkRJRFqVu/P80s385G8r2byrjIsmD+E754xlQI6mtxBpi5REpNV8sGU3s59bzvx1Oxg/OIdfXzaZKXm6h0ekLVMSkaTbVVrBf71cyP/M30CPLpncfuEELp2aS7qG44q0eUoikjRV1c4TCzdy54uF7Cw9wBUnDuffPjOanl21lKtIe6EkIkmxaMMnzH5uOe9v2sW0vN7MPn884wbnRB2WiLQwJRFpUcV7yvjZ3z/g6cWbGJiTxS8vzef8SYM17YhIO6UkIi3iQGU1D/3zQ3716hoOVFZz48yRfPW0o+nWWf/ERNoz/Q+XZnN3vvynd3h99TZOH9uf7583jhF9u0Udloi0AiURabb/W1TE66u38f3zxvHlU0ZEHY6ItCLdHizNsm1vObf/dSXT8nrzpZPzog5HRFqZkog0y38+v4L9B6r4yUUTNA27SAekJCIJKygsZs6Sj7nxtJEc3T876nBEJAJKIpKQ0gOVfO/ZZYzs142vzBwZdTgiEhF1rEtCfvHyKoo+2c8T159E5wyt+SHSUakmIk22bNMu/vDGh1w2LZdpIzSBokhHlnJJxMxmm9kmM1sSbv8Sc+wWM1tjZoVmdlaUcXZUlVXVfPfppfTp3pnvnjM26nBEJGKp2pz1C3e/K7bAzMYBlwLjgcHAK2Y22t2rogiwo3ron+tZtmk3933hOHp0yYw6HBGJWMrVRA5jFvCYu5e7+4fAGmBaxDF1KBt3lHL3S6v49Nj+/MuxA6MOR0RSQKomka+Z2VIze9DMeoVlQ4CNMecUhWWHMLPrzGyhmS0sKSlJdqwdgrvzvWeXYQY/umCCJlQUESCiJGJmr5jZsga2WcBvgZFAPrAZuLupr+/u97v7FHef0q9fvxaOvmP6y9LNzF1VwrfOHMOQnl2iDkdEUkQkfSLufkY855nZA8Dz4e4mYFjM4aFhmSTZztID/Ogvy5k0tAdXa2oTEYmRcs1ZZjYoZvdCYFn4/DngUjPrbGYjgFHAgtaOryP66d8+4JPSCn560UQtaSsiB4m7JmJmXd29NJnBhO4ws3zAgfXA9QDuvtzMngBWAJXAVzUyK/neWrudxxdu5IZTR2plQhE5xBGTiJmdDPwe6A7kmtkk4Hp3vzEZAbn7lYc5djtwezLeVw5VVlHFfzzzPrm9u3Lzp0dFHY6IpKB4mrN+AZwFbAdw9/eAGckMSlLDfa+t4cNt+7j9wgl06aSpTUTkUHH1ibj7xnpFakZq5wq37OG3BWu5aPIQpo/SCDcRaVg8fSIbwyYtN7NM4GZgZXLDkihVVzu3PL2U7KwMbj33mKjDEZEUFk9N5AbgqwQ39m0iuH8jKf0hkhoeXfARiz/ayffOHUef7p2jDkdEUlg8NZEx7n55bIGZfQp4MzkhSZS27Crjjr9/wKeO7sNFxzU4IYCISK14aiK/jrNM2oHZzy3nQFU1t19wrKY2EZEjarQmYmYnAScD/czs32IO5QAaqtMOvbh8Cy8s38K/nz2GvL7dog5HRNqAwzVndSK4NyQDiF1Aezfw2WQGJa1vT1kFt81ZztiB2Vw7/aiowxGRNqLRJOLuc4G5ZvaQu29oxZgkAne+WMjWPWX87srjyUxPudlwRCRFxdOxXmpmdxIsBpVVU+jupyctKmlVizZ8wv/M38DVJ+WRP6xn1OGISBsSz0/OR4EPgBHADwnms3oniTFJKzpQWc1/PP0+A3Oy+NZZY6IOR0TamHiSSB93/wNQ4e5z3f0aQLWQduKB19dRuHUP/zlrAt07p+pqySKSquL51qgIHzeb2bnAx0Dv5IUkreXDbfv45aur+ZdjB3LGuAFRhyMibVA8SeTHZtYD+CbB/SE5wDeSGpUklbvzyspifvK3lXTOSGP2v46POiQRaaOOmETcvWZlwV3AackNR5Jt4fod/OzvH7BwwyeM6NuN311xPP1zso58oYhIA+JZT2QE8HUgL/Z8dz8/eWFJS1u1dQ93vFDIKyu30i+7M7dfOIHPTRmm4bwi0izxNGc9C/wB+AtQndxwwMweB2qGCfUEdrp7vpnlEcweXBgem+/uNyQ7nrbu4537+cXLq3hqcRHdOmXwrTNHc80pI+jaSZ3oItJ88XyTlLn7r5IeScjdP1/z3MzuJmhGq7HW3fNbK5a2bGfpAX5TsJaH/rkeHK751AhuPO1oenfrFHVoItKOxJNEfmlmtwEvAeU1he6+OGlRARbM/vc5NJy4SfYfqOKP//yQ3xasZW95JRdNHso3PjOKob26Rh2aiLRD8SSRY4ErCb7Ma5qznOR/uU8Htrr76piyEWb2LsH8Xd9z99eTHEObUVlVzf8tKuKeV1axdXc5p4/tz7+fPYaxA3OiDk1E2rF4ksglwFHufqCl3tTMXgEGNnDoVnefEz6/DPhzzLHNQK67bzez44FnzWy8u+9u4PWvA64DyM3NbamwU5K78+Lyrdzx4gesK9nH5Nye/OrSyZxwVJ+oQxORDiCeJLKMoIO7uKXe1N3PONxxM8sALgKOj7mmnLA5zd0XmdlaYDSwsIHXvx+4H2DKlCneUnGnmvnrtvPzFz7g3Y92MrJfN/77yuM5c9wArQMiIq0mniTSE/jAzN7h4D6RZA7xPQP4wN2LagrMrB+ww92rzOwoYBSwLokxpKyVm3dzxwsf8FphCQNzsvj5xcdy8XFDydBwXRFpZfEkkduSHsWhLuXgpiyAGcCPzKyCoG/mBnff0eqRRezulwq597U1ZHfO4LvnjOWLJ+eRlak1wkQkGvHcsT63NQKp955fbKDsKeCp1o4llawp3st9r63hX44dxE8uOJYeXTOjDklEOrjDLY/7hrufYmZ7CEZj1R4C3N017KeV3fuP1XTOSOeH549XAhGRlHC4lQ1PCR+zGztHWs+a4r08997HXDv9KPp27xx1OCIiQBzriZjZ/8RTJslVUwu5dobWPxeR1BHPcJ6D5gkPh98e38i5kgRrS4JayFUnDVctRERSSqNJxMxuCftDJprZ7nDbA2wF5jR2nbS8e/+xRrUQEUlJjSYRd/9p2B9yp7vnhFu2u/dx91taMcYObW3JXuYs2aRaiIikpHias543s24AZnaFmf2XmQ1PclwSUi1ERFJZPEnkt0CpmU0iWCJ3LfBwUqMSoK4WcqVqISKSouJJIpXu7sAs4F53vw/QsN9WcO8/1tApI43rVAsRkRQVz7Qne8zsFoLp4KebWRqgO92SbF1YC/n/dF+IiKSweGoinyeYePEad98CDAXuTGpUolqIiLQJR0wiYeJ4Cqj5ObwNeCaZQXV060r28uySTVx1Up5qISKS0uK5Y/1a4Engv8OiIcCzyQyqo1MtRETainias74KfIpgSVrC5Wr7JzOojqymFnLliRqRJSKpL54kUh67NG447Um7XS0wanW1kJFRhyIickTxJJG5ZvYfQBcz+wzwf8BfkhtWxxRbC+mXrVqIiKS+eJLId4ES4H3geuBvwPeSGVRHde9rqoWISNsSz+isand/wN0vcffPhs+b3ZxlZpeY2XIzqzazKfWO3WJma8ys0MzOiik/OyxbY2bfbW4MqeTDbft49l3VQkSkbYmnJpIsy4CLgHmxhWY2jmCN9fHA2cBvzCzdzNKB+4BzgHHAZeG57cKv/7FatRARaXPiuWM9Kdx9JYCZ1T80C3jM3cuBD81sDTAtPLbG3deF1z0WnruidSJOnppayJdPGaFaiIi0KQnVRMzsrpYOJMYQYGPMflFY1lh5m6daiIi0VYk2Z30unpPM7BUzW9bANivB942LmV1nZgvNbGFJSUky36rZamohV5ygvhARaXsSbc46pA2qIe5+RgKvvQkYFrM/NCzjMOX13/d+4H6AKVOmpPQ9LbW1kFN1d7qItD2NJhEz693YIeJMIgl6DvhfM/svYDAwClgQvucoMxtBkDwuBb6QxDiSbv22fcxZ8jFfOjmP/tlZUYcjItJkh6uJLCK4M72hhFHR3Dc2swuBXwP9gL+a2RJ3P8vdl5vZEwQd5pXAV929Krzma8CLQDrwoLsvb24cUfr1P9aQmW6qhYhIm9VoEnH3Ecl8Y3d/hkZmA3b324HbGyj/G8HNjm3e+m37eHbJJtVCRKRNO1xz1nGHu9DdF7d8OB2HaiEi0h4crjlrIcENgdvC/dhmLQdOT1ZQ7V1NLeSLqoWISBt3uCTyb8Bngf3AY8Az7r63VaJq5379jzVkpBnXqxYiIm1co/eJuPs97n4K8HWCobWvmtkTZpbfatG1QzW1kCtOHK5aiIi0efFMwLgOmAO8RDD9yOhkB9We3fuaaiEi0n4crmP9KIJ7MWYRTDfyGPATd9/fSrG1O+u37eOZd9UXIiLtx+H6RNYASwlqIbuBXOArNRMmuvt/JT26dka1EBFpbw6XRH5E3TK43esdS+mpRFJRTS3k6pNUCxGR9uNwNxvObuyYmf2/pETTjtXUQm5QLURE2pFEZ/H9txaNop2rqYVcfsJw+ueoFiIi7UeiSSSZEzC2O/epFiIi7VSiSUR9InE6UFnN397fzAX5Q1QLEZF253BDfPfQcLIwoEvSImpnFm7Ywb4DVZx+TP+oQxERaXGH61jPbs1A2qu5hSVkphufOrpv1KGIiLS4RJuzJE5zV5UwZXhvundOdBFJEZHUpSSSRJt37eeDLXuYOaZf1KGIiCSFkkgSzS0sAWDmGPWHiEj7FEkSMbNLzGy5mVWb2ZSY8s+Y2SIzez98PD3mWIGZFZrZknBL+W/mgsISBuZkMXpA/Rv+RUTah6ga6pcBFwH/Xa98G/Cv7v6xmU0gWE99SMzxy919YSvF2CwVVdW8uWYb504cRM18YyIi7U0kScTdVwKHfLm6+7sxu8uBLmbW2d3LWzG8FrFowyfsKa9Uf4iItGup3CdyMbC4XgL5Y9iU9X07zM97M7vOzBaa2cKSkpLkR9qAuatKyEjT0F4Rad+SlkTM7BUzW9bANiuOa8cDPweujym+3N2PBaaH25WNXe/u97v7FHef0q9fNDWBgsISjh/ei+yszEjeX0SkNSStOcvdz0jkOjMbCjwDXOXua2Neb1P4uMfM/pdglcWHWyLWlrZ1dxkrN+/mO2ePjToUEZGkSqnmLDPrCfwV+K67vxlTnmFmfcPnmcB5BJ3zKalmaO+po9UfIiLtW1RDfC80syLgJOCvZvZieOhrwNHAD+oN5e0MvGhmS4ElwCbggShij0fBqmIG5HTmmEGaOUZE2reoRmc9Q9BkVb/8x8CPG7ns+KQG1UIqq6p5ffU2zpkwUEN7RaTdS6nmrPbg3Y072VNWqbvURaRDUBJpYQWFxaRraK+IdBBKIi2soLCE43N70aOLhvaKSPunJNKCiveUsfzj3Zyqu9RFpINQEmlBGtorIh2NkkgLKlhVQr/szowfnBN1KCIirUJJpIVUVlXzxuptnDq6n4b2ikiHoSTSQt4r2smu/RWatVdEOhQlkRZSUFhCmsH0o5VERKTjiGpRqnanoLCE43J70aOrhvaKtBUVFRUUFRVRVlYWdSiRycrKYujQoWRmJvbdpSTSAkr2lPP+pl188zOjow5FRJqgqKiI7Oxs8vLyOmRfpruzfft2ioqKGDFiREKvoeasFjBvVTC0V1OdiLQtZWVl9OnTp0MmEAhWl+3Tp0+zamJKIi1g7qoS+nbvpKG9Im1QR00gNZr7+ZVEmqmq2pm3uoQZo/uRltax/zGKSNN179496hCaRUmkmd4r2snO0go1ZYlISquqqkrK6yqJNFPN0N4ZozRrr4gkzt359re/zYQJEzj22GN5/PHHAfjqV7/Kc889B8CFF17INddcA8CDDz7IrbfeCsAjjzzCtGnTyM/P5/rrr69NGN27d+eb3/wmkyZN4q233kpK3JGMzjKzS4DZwDHANHdfGJbnASuBwvDU+e5+Q3jseOAhoAvwN+Bmd/fWjLshcwuLyR/Wk55dO0Udiog0ww//spwVH+9u0dccNziH2/51fFznPv300yxZsoT33nuPbdu2MXXqVGbMmMH06dN5/fXXOf/889m0aRObN28G4PXXX+fSSy9l5cqVPP7447z55ptkZmZy44038uijj3LVVVexb98+TjjhBO6+++4W/VyxoqqJLAMuAuY1cGytu+eH2w0x5b8FrgVGhdvZyQ/z8LbvLWfppl2cOlpNWSLSPG+88QaXXXYZ6enpDBgwgFNPPZV33nmnNomsWLGCcePGMWDAADZv3sxbb73FySefzKuvvsqiRYuYOnUq+fn5vPrqq6xbtw6A9PR0Lr744qTGHdXyuCsh/lEBZjYIyHH3+eH+w8AFwN+TFWM8Xl+9DXc01YlIOxBvjaG1DRkyhJ07d/LCCy8wY8YMduzYwRNPPEH37t3Jzs7G3bn66qv56U9/esi1WVlZpKenJzW+VOwTGWFm75rZXDObHpYNAYpizikKyyJVUFhMn26dOHZIj6hDEZE2bvr06Tz++ONUVVVRUlLCvHnzmDZtGgAnnngi99xzT23z1l133cX06cHX46c//WmefPJJiouLAdixYwcbNmxotbiTVhMxs1eAgQ0cutXd5zRy2WYg1923h30gz5pZk38emNl1wHUAubm5Tb08LtXVzrxw1l4N7RWR5rrwwgt56623mDRpEmbGHXfcwcCBwVfo9OnTeemllzj66KMZPnw4O3bsqE0i48aN48c//jFnnnkm1dXVZGaPRlKoAAAPmElEQVRmct999zF8+PBWidui7Js2swLgWzUd640dBzYBr7n72LD8MmCmu19/pPeYMmWKL1zY4Ms3y5KNO7ngvjf55aX5zMqPvFIkIglYuXIlxxxzTNRhRK6hPwczW+TuU450bUo1Z5lZPzNLD58fRdCBvs7dNwO7zexECzpSrgIaq820ioLCYsxg+ij1h4hIxxVJEjGzC82sCDgJ+KuZvRgemgEsNbMlwJPADe6+Izx2I/B7YA2wlog71QsKS5g4tCe9u2lor4h0XFGNznoGeKaB8qeApxq5ZiEwIcmhxeWTfQd4r2gnN50+KupQREQilVLNWW3FvNUlGtorIoKSSELmFpbQq2smE4f2jDoUEZFIKYk0UXW1M3dVMGtvuob2ikgHpyTSRMs+3sX2fQfUlCUigpJIkxUUlmAGMzS0V0TaGHenurq6RV9TSaSJ5q4q4dghPejTvXPUoYhIO7B+/XomTKgbeHrXXXcxe/ZsZs6cyc0330x+fj4TJkxgwYIFAMyePZsrr7ySk046iVGjRvHAAw/UXnvnnXcydepUJk6cyG233Vb7+mPGjOGqq65iwoQJbNy4sUXjj2SIb1u1s/QA7370CV877eioQxGRlvb378KW91v2NQceC+f8LOHLS0tLWbJkCfPmzeOaa65h2bJlACxdupT58+ezb98+Jk+ezLnnnsuyZctYvXo1CxYswN05//zzmTdvHrm5uaxevZo//elPnHjiiS31yWopiTTB66u3Ue1wqlYxFJFWcNlllwEwY8YMdu/ezc6dOwGYNWsWXbp0oUuXLpx22mksWLCAN954g5deeonJkycDsHfvXlavXk1ubi7Dhw9PSgIBJZEmKSgsoWfXTPKHaWivSLvTjBpDc2RkZBzUT1FWVlb7vP5yGTX7DZW7O7fccgvXX3/wlILr16+nW7duLR12LfWJxKlmaO/0URraKyItZ8CAARQXF7N9+3bKy8t5/vnna4/VLJH7xhtv0KNHD3r0CJadmDNnDmVlZWzfvp2CggKmTp3KWWedxYMPPsjevXsB2LRpU+308MmkmkicVmzezba95cwcrVFZItJyMjMz+cEPfsC0adMYMmQIY8eOrT2WlZXF5MmTqaio4MEHH6wtnzhxIqeddhrbtm3j+9//PoMHD2bw4MGsXLmSk046CQjWV3/kkUeSviiVkkic5q4qAWCGkoiItLCbbrqJm2666aCymTNncsUVV3DPPfcccv7EiRN5+OGHDym/+eabufnmmw8pr+mQTwY1Z8WpoLCYCUNy6Jetob0iIjVUE4nDrv0VLP5oJ185dWTUoYhIB1FQUNBg+ezZs1s1jiNRTSQOb6zeRlW1a6oTEZF6lETiUFBYTE5Whob2iojUoyRyBO7h0N7R/chI1x+XiEisqJbHvcTMlptZtZlNiSm/3MyWxGzVZpYfHisws8KYY61y2/jKzXso3qOhvSIiDYnqp/Uy4CJgXmyhuz/q7vnung9cCXzo7ktiTrm85ri7J/8uGqBgVfA2pyqJiEiS/OpXv+KYY47h8ssvb/Sc7t27A4dO2Bi1qNZYXwmH3rpfz2XAY60S0GEUFJYwblAO/XOyog5FRNqp3/zmN7zyyisMHTo06lCaLJUb+T8P/Lle2R/Dpqzv2xEyUEvYXVbBog2faFSWiCTNDTfcwLp16zjnnHPo0aMHd911V+2xCRMmsH79+uiCi0PSaiJm9gowsIFDt7r7nCNcewJQ6u6xt1le7u6bzCwbeIqguevQWzaD668DrgPIzc1NJHwA3qwd2qtZe0Xau1Vb97CnrKJFXzM7K5PRA7IPe87vfvc7XnjhBV577TXuvffeFn3/1pC0JOLuZzTj8kupVwtx903h4x4z+19gGo0kEXe/H7gfYMqUKZ5oEAWFJWRnZXBcrob2iog0JOXuWDezNOBzwPSYsgygp7tvM7NM4DzglWTGUTu0d1RfDe0V6QCOVGNoDYebFj5VRTXE90IzKwJOAv5qZi/GHJ4BbHT3dTFlnYEXzWwpsATYBDxAEhVu3cOW3WUalSUirSYvL4/FixcDsHjxYj788MOIIzqyqEZnPQM808ixAuDEemX7gOOTH1mdgsJg1t5TR6s/RERax8UXX8zDDz/M+PHjOeGEExg9enTUIR1RyjVnpYqCwmLGDsxmYA8N7RWR5IodgfXSSy81eE7NYlN5eXlJndq9qZREGuDujB/cg0FKICIih6Uk0gAz4/vnjYs6DBGRlKdhRyIikjAlERHp0NwTvpWsXWju51cSEZEOKysri+3bt3fYROLubN++naysxPt/1SciIh3W0KFDKSoqoqSkJOpQIpOVldWsiR+VRESkw8rMzGTEiBFRh9GmqTlLREQSpiQiIiIJUxIREZGEWXsflWBmJcCGqOMA+gLbog6iAYqraRRX0yiupkmluIa7+xFnoG33SSRVmNlCd58SdRz1Ka6mUVxNo7iaJlXjOhw1Z4mISMKUREREJGFKIq3n/qgDaITiahrF1TSKq2lSNa5GqU9EREQSppqIiIgkTEmkFZhZupm9a2bPRx1LLDPraWZPmtkHZrbSzE6KOiYAM/uGmS03s2Vm9mczi2R1MDN70MyKzWxZTFlvM3vZzFaHj71SJK47w7/HpWb2jJn1TIW4Yo5908zczPqmSlxm9vXwz2y5md2RCnGZWb6ZzTezJWa20MymtXZcTaUk0jpuBlZGHUQDfgm84O5jgUmkQIxmNgS4CZji7hOAdODSiMJ5CDi7Xtl3gVfdfRTwarjf2h7i0LheBia4+0RgFXBLawdFw3FhZsOAM4GPWjug0EPUi8vMTgNmAZPcfTxwVyrEBdwB/NDd84EfhPspTUkkycxsKHAu8PuoY4llZj2AGcAfANz9gLvvjDaqWhlAFzPLALoCH0cRhLvPA3bUK54F/Cl8/ifgglYNiobjcveX3L0y3J0PJD4tawvGFfoF8O9AJB2wjcT1FeBn7l4enlOcInE5kBM+70FE//abQkkk+e4h+A9UHXUg9YwASoA/hk1tvzezblEH5e6bCH4VfgRsBna5+0vRRnWQAe6+OXy+BRgQZTCNuAb4e9RBAJjZLGCTu78XdSz1jAamm9nbZjbXzKZGHVDo/wF3mtlGgv8HUdQom0RJJInM7Dyg2N0XRR1LAzKA44DfuvtkYB/RNM0cJOxjmEWQ5AYD3czsimijapgHQxtTanijmd0KVAKPpkAsXYH/IGiWSTUZQG/gRODbwBNmZtGGBAQ1pG+4+zDgG4QtBalMSSS5PgWcb2brgceA083skWhDqlUEFLn72+H+kwRJJWpnAB+6e4m7VwBPAydHHFOsrWY2CCB8bPVmkMaY2ReB84DLPTXG7o8k+DHwXvh/YCiw2MwGRhpVoAh42gMLCFoKWr3TvwFXE/ybB/g/QB3rHZm73+LuQ909j6Bz+B/unhK/qt19C7DRzMaERZ8GVkQYUo2PgBPNrGv4y/DTpECHf4znCP6jEz7OiTCWWmZ2NkGz6fnuXhp1PADu/r6793f3vPD/QBFwXPhvL2rPAqcBmNlooBOpMfHhx8Cp4fPTgdURxhIXrWzYsX0deNTMOgHrgC9FHA/u/raZPQksJmiWeZeI7uI1sz8DM4G+ZlYE3Ab8jKDp48sEs0N/LkXiugXoDLwctsrMd/cboo7L3SNvjmnkz+tB4MFweO0B4OrWrr01Ete1wC/DQSVlwHWtGVMidMe6iIgkTM1ZIiKSMCURERFJmJKIiIgkTElEREQSpiQiIiIJUxKRSIQzut4ds/8tM5vdQq/9kJl9tiVe6wjvc0k4+/Fr9coHh8OUm/JaXzSzexOM40dmdkYi1zZXPO9tZjPNLJVuGJUWpPtEJCrlwEVm9lN3T4WbvAAws4yYiQyP5MvAte7+Rmyhu38MJD2JxbxfZNOKxPneM4G9wD+TG41EQTURiUolwU2E36h/oH5Nwsz2ho8zw8ny5pjZOjP7mZldbmYLzOx9MxsZ8zJnhOsxrArnMKtZ1+VOM3snXHfj+pjXfd3MnqOBu/bN7LLw9ZeZ2c/Dsh8ApwB/MLM7652fV7NGRFjDeNrMXrBgDZI7Ys77UhjfAoIpcg77+cPn3wljec/Mflb/fDNbb2Y/NLPF4Xljw/J+Fqx/sjycbHODNbC2h5ntNbNfhOe9amb9wvKadS5q1ivpFc97m1kecAPwDQvWyJge1uCWhZ9hXv0YpG1REpEo3QdcbsG09PGaRPCldAxwJTDa3acRTLX/9Zjz8gjmHToX+J0FC1t9mWBW4KnAVOBaMxsRnn8ccLO7j459MzMbDPycYAqKfGCqmV3g7j8CFhLMU/XtI8ScD3weOBb4vJkNs2DerR8SJI9TgHFH+uBmdg7B5JQnuPskGl9rYpu7Hwf8FvhWWHYbwbQ74wnmSctt5NpuwMLwvLnhdQAPA98J1yt5P6b8sO/t7uuB3wG/cPd8d3+dYELGs8LPcP6RPrekNiURiYy77yb4crqpCZe94+6bw3Ug1gI108S/T5A4ajzh7tXuvppgSpexBAsjXWVmS4C3gT7AqPD8Be7+YQPvNxUoCCeErJkdd0YT4oVgEatd7l5GUNMZDpwQ87oHgMfjeJ0zgD/WzI3l7g2t3QF1E/gtou7P5BSCSUBx9xeATxq5tjomlkeAU8Ik39Pd54blf6LxP4OG3ru+N4GHzOxagkXHpA1TEpGo3UNQQ4hdy6SS8N+mmaURTI5XozzmeXXMfjUH9/HVn8/HAQO+Hv4iznf3ETFrlexr1qc4vNiYqzhyX+ThPn9T3i+e9zqSps6LdMT3Duf0+h4wDFhkZn0SD0+ipiQikQp/TT9BkEhqrAeOD5+fD2Qm8NKXmFla2E9yFFAIvAh8xcwyIZi91Y68ENcC4FQz62tm6cBlBM08zfV2+Lp9wnguiTm2noY//8vAlyxYpwMz692E93uTcLJIMzsTaGxt+DTqBgV8AXjD3XcBn5jZ9LD8Spr2Z7AHyK7ZMbOR7v522ClfQpBMpI3S6CxJBXcDX4vZfwCYY2bvAS+QWC3hI4IEkAPc4O5lZvZ7giaWxWZmBF9gh13e1t03m9l3gdcIajJ/dfdmT/8evu5s4C1gJ7Ak5nCDn9/dXzCzfGChmR0A/kaw6FM8fgj82cyuDN9zC8GXe337gGlm9j2CtVI+H5ZfTdC31JWmz/j8F+BJC1Y5/DpBJ/sogj/PV4FUW/VQmkCz+Ip0AGbWGahy90ozO4lgRcv8Bs7b6+7dWz9CaatUExHpGHIJ1kFJI1g/49qI45F2QjURERFJmDrWRUQkYUoiIiKSMCURERFJmJKIiIgkTElEREQSpiQiIiIJ+/8BTwaCL/b58K4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc35759cbe0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(fMs, fvfe_lml, label=\"lower\")\n",
    "plt.plot(fMs, fvupper_lml, label=\"upper\")\n",
    "plt.axhline(full_lml, label=\"full\", alpha=0.3)\n",
    "plt.xlabel(\"Number of inducing points\")\n",
    "plt.ylabel(\"LML estimate\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Now, as the hyperparameters are fixed, the bound _does_ monotonically decrease. We chose the optimal hyperparameters here, but the picture should be the same for any hyperparameter setting. This shows that we increasingly get a better estimate of the marginal likelihood as we add more inducing points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### A tight estimate bound does not imply a converged model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T12:33:03.794371Z",
     "start_time": "2018-06-19T12:32:58.518944Z"
    },
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lower bound: -62.487090\n",
      "Upper bound: -62.484641\n"
     ]
    }
   ],
   "source": [
    "vfe = gpflow.models.SGPR(X, Y, gpflow.kernels.RBF(1), X[None, 0, :].copy())\n",
    "vfe.compile()\n",
    "gpflow.train.ScipyOptimizer().minimize(vfe, maxiter=notebook_niter(1000))\n",
    "print(\"Lower bound: %f\" % vfe.compute_log_likelihood())\n",
    "print(\"Upper bound: %f\" % vfe.compute_upper_bound())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "In this case we show that for the hyperparameter setting, the bound is very tight. However, this does _not_ imply that we have enough inducing points, but simply that we have correctly identified the marginal likelihood for this particular hyperparameter setting. In this specific case, where we used a single inducing point, the model collapses to not using the GP at all (lengthscale is really long to model only the mean). The rest of the variance is explained by noise. This GP can be perfectly approximated with a single inducing point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T12:33:03.819380Z",
     "start_time": "2018-06-19T12:33:03.795936Z"
    },
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>class</th>\n",
       "      <th>prior</th>\n",
       "      <th>transform</th>\n",
       "      <th>trainable</th>\n",
       "      <th>shape</th>\n",
       "      <th>fixed_shape</th>\n",
       "      <th>value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>SGPR/likelihood/variance</th>\n",
       "      <td>Parameter</td>\n",
       "      <td>None</td>\n",
       "      <td>+ve</td>\n",
       "      <td>True</td>\n",
       "      <td>()</td>\n",
       "      <td>True</td>\n",
       "      <td>0.6824319445027769</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SGPR/feature/Z</th>\n",
       "      <td>Parameter</td>\n",
       "      <td>None</td>\n",
       "      <td>(none)</td>\n",
       "      <td>True</td>\n",
       "      <td>(1, 1)</td>\n",
       "      <td>True</td>\n",
       "      <td>[[2.6574094508]]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SGPR/kern/lengthscales</th>\n",
       "      <td>Parameter</td>\n",
       "      <td>None</td>\n",
       "      <td>+ve</td>\n",
       "      <td>True</td>\n",
       "      <td>()</td>\n",
       "      <td>True</td>\n",
       "      <td>1014.3847523740365</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SGPR/kern/variance</th>\n",
       "      <td>Parameter</td>\n",
       "      <td>None</td>\n",
       "      <td>+ve</td>\n",
       "      <td>True</td>\n",
       "      <td>()</td>\n",
       "      <td>True</td>\n",
       "      <td>0.10775256810399733</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                              class prior transform  trainable   shape  \\\n",
       "SGPR/likelihood/variance  Parameter  None       +ve       True      ()   \n",
       "SGPR/feature/Z            Parameter  None    (none)       True  (1, 1)   \n",
       "SGPR/kern/lengthscales    Parameter  None       +ve       True      ()   \n",
       "SGPR/kern/variance        Parameter  None       +ve       True      ()   \n",
       "\n",
       "                          fixed_shape                value  \n",
       "SGPR/likelihood/variance         True   0.6824319445027769  \n",
       "SGPR/feature/Z                   True     [[2.6574094508]]  \n",
       "SGPR/kern/lengthscales           True   1014.3847523740365  \n",
       "SGPR/kern/variance               True  0.10775256810399733  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vfe.as_pandas_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "This can be diagnosed by showing that there are other hyperparameter settings with higher upper bounds. This indicates that there may be better hyperparameter settings, but we cannot identify them due to the lack of inducing points. An example of this can be seen in the previous section."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
